MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
imsreorderingmodule Module Reference

Functions/Subroutines

subroutine, public ims_odrv (n, nja, nsp, ia, ja, p, ip, isp, flag)
 
subroutine ims_md (n, nja, ia, ja, mmax, v, l, head, last, next, mark, flag)
 
subroutine ims_mdi (n, nja, ia, ja, mmax, v, l, head, last, next, mark, tag, flag)
 
subroutine ims_mdm (n, mmax, vk, tail, v, l, last, next, mark)
 
subroutine ims_mdp (n, mmax, k, ek, tail, v, l, head, last, next, mark)
 
subroutine ims_mdu (n, mmax, ek, dmin, v, l, head, last, next, mark)
 

Function/Subroutine Documentation

◆ ims_md()

subroutine imsreorderingmodule::ims_md ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  nja,
integer(i4b), dimension(n + 1), intent(in)  ia,
integer(i4b), dimension(nja), intent(in)  ja,
integer(i4b), intent(in)  mmax,
integer(i4b), dimension(mmax), intent(inout)  v,
integer(i4b), dimension(mmax), intent(inout)  l,
integer(i4b), dimension(n), intent(inout)  head,
integer(i4b), dimension(n), intent(inout)  last,
integer(i4b), dimension(n), intent(inout)  next,
integer(i4b), dimension(n), intent(inout)  mark,
integer(i4b), intent(inout)  flag 
)
private

Definition at line 198 of file ImsReordering.f90.

200  !
201  !*****************************************************************
202  ! ims_md -- minimum degree algorithm (based on element model)
203  !*****************************************************************
204  !
205  ! description
206  !
207  ! ims_md finds a minimum degree ordering of the rows and
208  ! columns of a general sparse matrix m stored in (ia,ja,a)
209  ! format. when the structure of m is nonsymmetric, the ordering
210  ! is that obtained for the symmetric matrix m + m-transpose.
211  !
212  !
213  ! additional parameters
214  !
215  ! mmax - declared dimension of the one-dimensional arrays v and l.
216  ! mmax must be at least n+2k, where k is the number of
217  ! nonzeroes in the strict upper triangle of m
218  !
219  ! v - integer one-dimensional work array. dimension = mmax
220  !
221  ! l - integer one-dimensional work array. dimension = mmax
222  !
223  ! head - integer one-dimensional work array. dimension = n
224  !
225  ! last - integer one-dimensional array used to return the permutation
226  ! of the rows and columns of m corresponding to the minimum
227  ! degree ordering. dimension = n
228  !
229  ! next - integer one-dimensional array used to return the inverse of
230  ! the permutation returned in last. dimension = n
231  !
232  ! mark - integer one-dimensional work array (may be the same as v).
233  ! dimension = n
234  !
235  ! flag - integer error flag. values and their meanings are -
236  ! 0 no errors detected
237  ! 11n+1 insufficient storage in md
238  !
239  !
240  ! definitions of internal parameters
241  !
242  ! ---------+---------------------------------------------------------
243  ! v(s) - value field of list entry
244  ! ---------+---------------------------------------------------------
245  ! l(s) - link field of list entry (0 =) end of list)
246  ! ---------+---------------------------------------------------------
247  ! l(vi) - pointer to element list of uneliminated vertex vi
248  ! ---------+---------------------------------------------------------
249  ! l(ej) - pointer to boundary list of active element ej
250  ! ---------+---------------------------------------------------------
251  ! head(d) - vj =) vj head of d-list d
252  ! - 0 =) no vertex in d-list d
253  !
254  !
255  ! - vi uneliminated vertex
256  ! - vi in ek - vi not in ek
257  ! ---------+-----------------------------+---------------------------
258  ! next(vi) - undefined but nonnegative - vj =) vj next in d-list
259  ! - - 0 =) vi tail of d-list
260  ! ---------+-----------------------------+---------------------------
261  ! last(vi) - (not set until mdp) - -d =) vi head of d-list d
262  ! --vk =) compute degree - vj =) vj last in d-list
263  ! - ej =) vi prototype of ej - 0 =) vi not in any d-list
264  ! - 0 =) do not compute degree -
265  ! ---------+-----------------------------+---------------------------
266  ! mark(vi) - mark(vk) - nonneg. tag .lt. mark(vk)
267  !
268  !
269  ! - vi eliminated vertex
270  ! - ei active element - otherwise
271  ! ---------+-----------------------------+---------------------------
272  ! next(vi) - -j =) vi was j-th vertex - -j =) vi was j-th vertex
273  ! - to be eliminated - to be eliminated
274  ! ---------+-----------------------------+---------------------------
275  ! last(vi) - m =) size of ei = m - undefined
276  ! ---------+-----------------------------+---------------------------
277  ! mark(vi) - -m =) overlap count of ei - undefined
278  ! - with ek = m -
279  ! - otherwise nonnegative tag -
280  ! - .lt. mark(vk) -
281  !
282  !-----------------------------------------------------------------------
283  !
284  implicit none
285 
286  ! -- dummy variables
287  integer(I4B), intent(in) :: n
288  integer(I4B), intent(in) :: nja
289  integer(I4B), dimension(n + 1), intent(in) :: ia
290  integer(I4B), dimension(nja), intent(in) :: ja
291  integer(I4B), intent(in) :: mmax
292  integer(I4B), dimension(mmax), intent(inout) :: v
293  integer(I4B), dimension(mmax), intent(inout) :: l
294  integer(I4B), dimension(n), intent(inout) :: head
295  integer(I4B), dimension(n), intent(inout) :: last
296  integer(I4B), dimension(n), intent(inout) :: next
297  integer(I4B), dimension(n), intent(inout) :: mark
298  integer(I4B), intent(inout) :: flag
299 
300  ! -- local
301  integer(I4B) :: tag
302  integer(I4B) :: dmin
303  integer(I4B), pointer :: vk
304  integer(I4B), pointer :: ek
305  integer(I4B) :: tail
306  integer(I4B) :: k
307 
308  allocate (vk)
309  ek => vk
310  !
311  ! initialization
312  tag = 0
313  call ims_mdi(n, nja, ia, ja, mmax, v, l, head, last, next, &
314  mark, tag, flag)
315  if (flag .ne. 0) return
316  !
317  k = 0
318  dmin = 1
319  !
320  ! while k .lt. n do
321 1 if (k >= n) go to 4
322  !
323  ! search for vertex of minimum degree
324 2 if (head(dmin) > 0) go to 3
325  dmin = dmin + 1
326  go to 2
327  !
328  ! remove vertex vk of minimum degree from degree list
329 3 vk = head(dmin)
330  head(dmin) = next(vk)
331  if (head(dmin) > 0) last(head(dmin)) = -dmin
332  !
333  ! number vertex vk, adjust tag, and tag vk
334  k = k + 1
335  next(vk) = -k
336  last(ek) = dmin - 1
337  tag = tag + last(ek)
338  mark(vk) = tag
339  !
340  ! form element ek from uneliminated neighbors of vk
341  call ims_mdm(n, mmax, vk, tail, v, l, last, next, mark)
342  !
343  ! purge inactive elements and do mass elimination
344  call ims_mdp(n, mmax, k, ek, tail, v, l, head, last, next, mark)
345  !
346  ! update degrees of uneliminated vertices in ek
347  call ims_mdu(n, mmax, ek, dmin, v, l, head, last, next, mark)
348  !
349  go to 1
350  !
351  ! generate inverse permutation from permutation
352 4 do k = 1, n
353  next(k) = -next(k)
354  last(next(k)) = k
355  end do
356  !
357  deallocate (vk)
358  !
359  return
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_mdi()

subroutine imsreorderingmodule::ims_mdi ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  nja,
integer(i4b), dimension(n + 1), intent(in)  ia,
integer(i4b), dimension(nja), intent(in)  ja,
integer(i4b), intent(in)  mmax,
integer(i4b), dimension(mmax), intent(inout)  v,
integer(i4b), dimension(mmax), intent(inout)  l,
integer(i4b), dimension(n), intent(inout)  head,
integer(i4b), dimension(n), intent(inout)  last,
integer(i4b), dimension(n), intent(inout)  next,
integer(i4b), dimension(n), intent(inout)  mark,
integer(i4b), intent(in)  tag,
integer(i4b), intent(inout)  flag 
)
private

Definition at line 362 of file ImsReordering.f90.

364  !
365  !***********************************************************************
366  ! ims_mdi -- initialization
367  !***********************************************************************
368  implicit none
369 
370  ! -- dummy variables
371  integer(I4B), intent(in) :: n
372  integer(I4B), intent(in) :: nja
373  integer(I4B), dimension(n + 1), intent(in) :: ia
374  integer(I4B), dimension(nja), intent(in) :: ja
375  integer(I4B), intent(in) :: mmax
376  integer(I4B), dimension(mmax), intent(inout) :: v
377  integer(I4B), dimension(mmax), intent(inout) :: l
378  integer(I4B), dimension(n), intent(inout) :: head
379  integer(I4B), dimension(n), intent(inout) :: last
380  integer(I4B), dimension(n), intent(inout) :: next
381  integer(I4B), dimension(n), intent(inout) :: mark
382  integer(I4B), intent(in) :: tag
383  integer(I4B), intent(inout) :: flag
384 
385  ! -- local
386  integer(I4B) :: sfs
387  integer(I4B) :: vi
388  integer(I4B) :: dvi
389  integer(I4B) :: vj
390  integer(I4B) :: jmin
391  integer(I4B) :: jmax
392  integer(I4B) :: j
393  integer(I4B) :: lvk
394  integer(I4B) :: kmax
395  integer(I4B) :: k
396  integer(I4B) :: nextvi
397  integer(I4B) :: ieval
398  !
399  ! initialize degrees, element lists, and degree lists
400  do vi = 1, n
401  mark(vi) = 1
402  l(vi) = 0
403  head(vi) = 0
404  end do
405  sfs = n + 1
406  !
407  ! create nonzero structure
408  ! for each nonzero entry a(vi,vj)
409  louter: do vi = 1, n
410  jmin = ia(vi)
411  jmax = ia(vi + 1) - 1
412  if (jmin > jmax) cycle louter
413  linner1: do j = jmin, jmax !5
414  vj = ja(j)
415  !if (vj-vi) 2, 5, 4
416  ieval = vj - vi
417  if (ieval == 0) cycle linner1 !5
418  if (ieval > 0) go to 4
419  !
420  ! if a(vi,vj) is in strict lower triangle
421  ! check for previous occurrence of a(vj,vi)
422  lvk = vi
423  kmax = mark(vi) - 1
424  if (kmax == 0) go to 4
425  linner2: do k = 1, kmax
426  lvk = l(lvk)
427  if (v(lvk) == vj) cycle linner1 !5
428  end do linner2
429  ! for unentered entries a(vi,vj)
430 4 if (sfs >= mmax) go to 101
431  !
432  ! enter vj in element list for vi
433  mark(vi) = mark(vi) + 1
434  v(sfs) = vj
435  l(sfs) = l(vi)
436  l(vi) = sfs
437  sfs = sfs + 1
438  !
439  ! enter vi in element list for vj
440  mark(vj) = mark(vj) + 1
441  v(sfs) = vi
442  l(sfs) = l(vj)
443  l(vj) = sfs
444  sfs = sfs + 1
445  end do linner1
446  end do louter
447  !
448  ! create degree lists and initialize mark vector
449  do vi = 1, n
450  dvi = mark(vi)
451  next(vi) = head(dvi)
452  head(dvi) = vi
453  last(vi) = -dvi
454  nextvi = next(vi)
455  if (nextvi > 0) last(nextvi) = vi
456  mark(vi) = tag
457  end do
458  !
459  return
460  !
461  ! ** error- insufficient storage
462 101 flag = 9 * n + vi
463  return
Here is the caller graph for this function:

◆ ims_mdm()

subroutine imsreorderingmodule::ims_mdm ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  mmax,
integer(i4b), intent(in)  vk,
integer(i4b), intent(inout)  tail,
integer(i4b), dimension(mmax), intent(inout)  v,
integer(i4b), dimension(mmax), intent(inout)  l,
integer(i4b), dimension(n), intent(inout)  last,
integer(i4b), dimension(n), intent(inout)  next,
integer(i4b), dimension(n), intent(inout)  mark 
)
private

Definition at line 466 of file ImsReordering.f90.

467  !
468  !***********************************************************************
469  ! ims_mdm -- form element from uneliminated neighbors of vk
470  !***********************************************************************
471  implicit none
472 
473  ! -- dummy variables
474  integer(I4B), intent(in) :: n
475  integer(I4B), intent(in) :: mmax
476  integer(I4B), intent(in) :: vk
477  integer(I4B), intent(inout) :: tail
478  integer(I4B), dimension(mmax), intent(inout) :: v
479  integer(I4B), dimension(mmax), intent(inout) :: l
480  integer(I4B), dimension(n), intent(inout) :: last
481  integer(I4B), dimension(n), intent(inout) :: next
482  integer(I4B), dimension(n), intent(inout) :: mark
483 
484  ! -- local
485  integer(I4B) :: tag
486  integer(I4B) :: s
487  integer(I4B) :: ls
488  integer(I4B), pointer :: vs
489  integer(I4B), pointer :: es
490  integer(I4B) :: b
491  integer(I4B) :: lb
492  integer(I4B) :: vb
493  integer(I4B) :: blp
494  integer(I4B) :: blpmax
495 
496  allocate (vs)
497  es => vs
498  !
499  ! initialize tag and list of uneliminated neighbors
500  tag = mark(vk)
501  tail = vk
502  !
503  ! for each vertex/element vs/es in element list of vk
504  ls = l(vk)
505 1 s = ls
506  if (s == 0) go to 5
507  ls = l(s)
508  vs = v(s)
509  if (next(vs) < 0) go to 2
510  !
511  ! if vs is uneliminated vertex, then tag and append to list of
512  ! uneliminated neighbors
513  mark(vs) = tag
514  l(tail) = s
515  tail = s
516  go to 4
517  !
518  ! if es is active element, then ...
519  ! for each vertex vb in boundary list of element es
520 2 lb = l(es)
521  blpmax = last(es)
522  louter: do blp = 1, blpmax !3
523  b = lb
524  lb = l(b)
525  vb = v(b)
526  !
527  ! if vb is untagged vertex, then tag and append to list of
528  ! uneliminated neighbors
529  if (mark(vb) >= tag) cycle louter !3
530  mark(vb) = tag
531  l(tail) = b
532  tail = b
533  end do louter
534  !
535  ! mark es inactive
536  mark(es) = tag
537  !
538 4 go to 1
539  !
540  ! terminate list of uneliminated neighbors
541 5 l(tail) = 0
542  !
543  deallocate (vs)
544  !
545  return
Here is the caller graph for this function:

◆ ims_mdp()

subroutine imsreorderingmodule::ims_mdp ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  mmax,
integer(i4b), intent(inout)  k,
integer(i4b), intent(in)  ek,
integer(i4b), intent(inout)  tail,
integer(i4b), dimension(mmax), intent(inout)  v,
integer(i4b), dimension(mmax), intent(inout)  l,
integer(i4b), dimension(n), intent(inout)  head,
integer(i4b), dimension(n), intent(inout)  last,
integer(i4b), dimension(n), intent(inout)  next,
integer(i4b), dimension(n), intent(inout)  mark 
)
private

Definition at line 548 of file ImsReordering.f90.

549  !
550  !***********************************************************************
551  ! ims_mdp -- purge inactive elements and do mass elimination
552  !***********************************************************************
553  implicit none
554 
555  ! -- dummy variables
556  integer(I4B), intent(in) :: n
557  integer(I4B), intent(in) :: mmax
558  integer(I4B), intent(inout) :: k
559  integer(I4B), intent(in) :: ek
560  integer(I4B), intent(inout) :: tail
561  integer(I4B), dimension(mmax), intent(inout) :: v
562  integer(I4B), dimension(mmax), intent(inout) :: l
563  integer(I4B), dimension(n), intent(inout) :: head
564  integer(I4B), dimension(n), intent(inout) :: last
565  integer(I4B), dimension(n), intent(inout) :: next
566  integer(I4B), dimension(n), intent(inout) :: mark
567 
568  ! -- local
569  integer(I4B) :: tag
570  integer(I4B) :: free
571  integer(I4B) :: li
572  integer(I4B) :: vi
573  integer(I4B) :: lvi
574  integer(I4B) :: evi
575  integer(I4B) :: s
576  integer(I4B) :: ls
577  integer(I4B) :: es
578  integer(I4B) :: ilp
579  integer(I4B) :: ilpmax
580  integer(I4B) :: i
581  !
582  ! initialize tag
583  tag = mark(ek)
584  !
585  ! for each vertex vi in ek
586  li = ek
587  ilpmax = last(ek)
588  if (ilpmax <= 0) go to 12
589  louter: do ilp = 1, ilpmax !11
590  i = li
591  li = l(i)
592  vi = v(li)
593  !
594  ! remove vi from degree list
595  if (last(vi) == 0) go to 3
596  if (last(vi) > 0) go to 1
597  head(-last(vi)) = next(vi)
598  go to 2
599 1 next(last(vi)) = next(vi)
600 2 if (next(vi) > 0) last(next(vi)) = last(vi)
601  !
602  ! remove inactive items from element list of vi
603 3 ls = vi
604 4 s = ls
605  ls = l(s)
606  if (ls == 0) go to 6
607  es = v(ls)
608  if (mark(es) < tag) go to 5
609  free = ls
610  l(s) = l(ls)
611  ls = s
612 5 go to 4
613  !
614  ! if vi is interior vertex, then remove from list and eliminate
615 
616 6 lvi = l(vi)
617  if (lvi .ne. 0) go to 7
618  l(i) = l(li)
619  li = i
620  !
621  k = k + 1
622  next(vi) = -k
623  last(ek) = last(ek) - 1
624  cycle louter !11
625  !
626  ! else ...
627  ! classify vertex vi
628 7 if (l(lvi) .ne. 0) go to 9
629  evi = v(lvi)
630  if (next(evi) >= 0) go to 9
631  if (mark(evi) < 0) go to 8
632  !
633  ! if vi is prototype vertex, then mark as such, initialize
634  ! overlap count for corresponding element, and move vi to end
635  ! of boundary list
636  last(vi) = evi
637  mark(evi) = -1
638  l(tail) = li
639  tail = li
640  l(i) = l(li)
641  li = i
642  go to 10
643  !
644  ! else if vi is duplicate vertex, then mark as such and adjust
645  ! overlap count for corresponding element
646 8 last(vi) = 0
647  mark(evi) = mark(evi) - 1
648  go to 10
649  !
650  ! else mark vi to compute degree
651 9 last(vi) = -ek
652  !
653  ! insert ek in element list of vi
654 10 v(free) = ek
655  l(free) = l(vi)
656  l(vi) = free
657  end do louter !11
658  !
659  ! terminate boundary list
660 12 l(tail) = 0
661  !
662  return
Here is the caller graph for this function:

◆ ims_mdu()

subroutine imsreorderingmodule::ims_mdu ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  mmax,
integer(i4b), intent(in)  ek,
integer(i4b), intent(inout)  dmin,
integer(i4b), dimension(mmax), intent(inout)  v,
integer(i4b), dimension(mmax), intent(inout)  l,
integer(i4b), dimension(n), intent(inout)  head,
integer(i4b), dimension(n), intent(inout)  last,
integer(i4b), dimension(n), intent(inout)  next,
integer(i4b), dimension(n), intent(inout)  mark 
)
private

Definition at line 665 of file ImsReordering.f90.

666  !
667  !***********************************************************************
668  ! ims_mdu -- update degrees of uneliminated vertices in ek
669  !***********************************************************************
670  implicit none
671 
672  ! -- dummy variables
673  integer(I4B), intent(in) :: n
674  integer(I4B), intent(in) :: mmax
675  integer(I4B), intent(in) :: ek
676  integer(I4B), intent(inout) :: dmin
677  integer(I4B), dimension(mmax), intent(inout) :: v
678  integer(I4B), dimension(mmax), intent(inout) :: l
679  integer(I4B), dimension(n), intent(inout) :: head
680  integer(I4B), dimension(n), intent(inout) :: last
681  integer(I4B), dimension(n), intent(inout) :: next
682  integer(I4B), dimension(n), intent(inout) :: mark
683 
684  ! -- local
685  integer(I4B) :: tag
686  integer(I4B) :: vi
687  integer(I4B) :: evi
688  integer(I4B) :: dvi
689  integer(I4B) :: s
690  integer(I4B), pointer :: vs
691  integer(I4B), pointer :: es
692  integer(I4B) :: b
693  integer(I4B) :: vb
694  integer(I4B) :: ilp
695  integer(I4B) :: ilpmax
696  integer(I4B) :: blp
697  integer(I4B) :: blpmax
698  integer(I4B) :: i
699 
700  allocate (vs)
701  es => vs
702  !
703  ! initialize tag
704  tag = mark(ek) - last(ek)
705  !
706  ! for each vertex vi in ek
707  i = ek
708  ilpmax = last(ek)
709  if (ilpmax <= 0) go to 11
710  louter: do ilp = 1, ilpmax !10
711  i = l(i)
712  vi = v(i)
713  !if (last(vi)) 1, 10, 8
714  if (last(vi) == 0) cycle louter !10
715  if (last(vi) > 0) goto 8
716  !
717  ! if vi neither prototype nor duplicate vertex, then merge elements
718  ! to compute degree
719  tag = tag + 1
720  dvi = last(ek)
721  !
722  ! for each vertex/element vs/es in element list of vi
723  s = l(vi)
724 2 s = l(s)
725  if (s == 0) go to 9
726  vs = v(s)
727  if (next(vs) < 0) go to 3
728  !
729  ! if vs is uneliminated vertex, then tag and adjust degree
730  mark(vs) = tag
731  dvi = dvi + 1
732  go to 5
733  !
734  ! if es is active element, then expand
735  ! check for outmatched vertex
736 3 if (mark(es) < 0) go to 6
737  !
738  ! for each vertex vb in es
739  b = es
740  blpmax = last(es)
741  linner: do blp = 1, blpmax !4
742  b = l(b)
743  vb = v(b)
744  !
745  ! if vb is untagged, then tag and adjust degree
746  if (mark(vb) >= tag) cycle linner !4
747  mark(vb) = tag
748  dvi = dvi + 1
749  end do linner !4
750  !
751 5 go to 2
752  !
753  ! else if vi is outmatched vertex, then adjust overlaps but do not
754  ! compute degree
755 6 last(vi) = 0
756  mark(es) = mark(es) - 1
757 7 s = l(s)
758  if (s == 0) cycle louter !10
759  es = v(s)
760  if (mark(es) < 0) mark(es) = mark(es) - 1
761  go to 7
762  !
763  ! else if vi is prototype vertex, then calculate degree by
764  ! inclusion/exclusion and reset overlap count
765 8 evi = last(vi)
766  dvi = last(ek) + last(evi) + mark(evi)
767  mark(evi) = 0
768  !
769  ! insert vi in appropriate degree list
770 9 next(vi) = head(dvi)
771  head(dvi) = vi
772  last(vi) = -dvi
773  if (next(vi) > 0) last(next(vi)) = vi
774  if (dvi < dmin) dmin = dvi
775  !
776  end do louter !10
777  !
778 11 deallocate (vs)
779  !
780  return
Here is the caller graph for this function:

◆ ims_odrv()

subroutine, public imsreorderingmodule::ims_odrv ( integer(i4b), intent(in)  n,
integer(i4b), intent(in)  nja,
integer(i4b), intent(in)  nsp,
integer(i4b), dimension(n + 1), intent(in)  ia,
integer(i4b), dimension(nja), intent(in)  ja,
integer(i4b), dimension(n), intent(inout)  p,
integer(i4b), dimension(n), intent(inout)  ip,
integer(i4b), dimension(nsp), intent(inout)  isp,
integer(i4b), intent(inout)  flag 
)

Definition at line 7 of file ImsReordering.f90.

8  !
9  ! 3/12/82
10  !***********************************************************************
11  ! odrv -- driver for sparse matrix reordering routines
12  !***********************************************************************
13  !
14  ! description
15  !
16  ! odrv finds a minimum degree ordering of the rows and columns
17  ! of a matrix m stored in (ia,ja,a) format (see below). for the
18  ! reordered matrix, the work and storage required to perform
19  ! gaussian elimination is (usually) significantly less.
20  !
21  ! note.. odrv and its subordinate routines have been modified to
22  ! compute orderings for general matrices, not necessarily having any
23  ! symmetry. the minimum degree ordering is computed for the
24  ! structure of the symmetric matrix m + m-transpose.
25  ! modifications to the original odrv module have been made in
26  ! the coding in subroutine mdi, and in the initial comments in
27  ! subroutines odrv and md.
28  !
29  ! if only the nonzero entries in the upper triangle of m are being
30  ! stored, then odrv symmetrically reorders (ia,ja,a), (optionally)
31  ! with the diagonal entries placed first in each row. this is to
32  ! ensure that if m(i,j) will be in the upper triangle of m with
33  ! respect to the new ordering, then m(i,j) is stored in row i (and
34  ! thus m(j,i) is not stored), whereas if m(i,j) will be in the
35  ! strict lower triangle of m, then m(j,i) is stored in row j (and
36  ! thus m(i,j) is not stored).
37  !
38  !
39  ! storage of sparse matrices
40  !
41  ! the nonzero entries of the matrix m are stored row-by-row in the
42  ! array a. to identify the individual nonzero entries in each row,
43  ! we need to know in which column each entry lies. these column
44  ! indices are stored in the array ja. i.e., if a(k) = m(i,j), then
45  ! ja(k) = j. to identify the individual rows, we need to know where
46  ! each row starts. these row pointers are stored in the array ia.
47  ! i.e., if m(i,j) is the first nonzero entry (stored) in the i-th row
48  ! and a(k) = m(i,j), then ia(i) = k. moreover, ia(n+1) points to
49  ! the first location following the last element in the last row.
50  ! thus, the number of entries in the i-th row is ia(i+1) - ia(i),
51  ! the nonzero entries in the i-th row are stored consecutively in
52  !
53  ! a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
54  !
55  ! and the corresponding column indices are stored consecutively in
56  !
57  ! ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
58  !
59  ! since the coefficient matrix is symmetric, only the nonzero entries
60  ! in the upper triangle need be stored. for example, the matrix
61  !
62  ! ( 1 0 2 3 0 )
63  ! ( 0 4 0 0 0 )
64  ! m = ( 2 0 5 6 0 )
65  ! ( 3 0 6 7 8 )
66  ! ( 0 0 0 8 9 )
67  !
68  ! could be stored as
69  !
70  ! - 1 2 3 4 5 6 7 8 9 10 11 12 13
71  ! ---+--------------------------------------
72  ! ia - 1 4 5 8 12 14
73  ! ja - 1 3 4 2 1 3 4 1 3 4 5 4 5
74  ! a - 1 2 3 4 2 5 6 3 6 7 8 8 9
75  !
76  ! or (symmetrically) as
77  !
78  ! - 1 2 3 4 5 6 7 8 9
79  ! ---+--------------------------
80  ! ia - 1 4 5 7 9 10
81  ! ja - 1 3 4 2 3 4 4 5 5
82  ! a - 1 2 3 4 5 6 7 8 9 .
83  !
84  !
85  ! parameters
86  !
87  ! n - order of the matrix
88  !
89  ! nja - number of nonzeroes in the matrix
90  !
91  ! nsp - declared dimension of the one-dimensional array isp. nsp
92  ! must be at least 3n+4k, where k is the number of nonzeroes
93  ! in the strict upper triangle of m
94  !
95  ! ia - integer one-dimensional array containing pointers to delimit
96  ! rows in ja and a. dimension = n+1
97  !
98  ! ja - integer one-dimensional array containing the column indices
99  ! corresponding to the elements of a. dimension = number of
100  ! nonzero entries in (the upper triangle of) m
101  !
102  ! a - real one-dimensional array containing the nonzero entries in
103  ! (the upper triangle of) m, stored by rows. dimension =
104  ! number of nonzero entries in (the upper triangle of) m
105  !
106  ! p - integer one-dimensional array used to return the permutation
107  ! of the rows and columns of m corresponding to the minimum
108  ! degree ordering. dimension = n
109  !
110  ! ip - integer one-dimensional array used to return the inverse of
111  ! the permutation returned in p. dimension = n
112  !
113  ! isp - integer one-dimensional array used for working storage.
114  ! dimension = nsp
115  !
116  ! path - integer path specification. values and their meanings are -
117  ! 1 find minimum degree ordering only
118  ! 2 find minimum degree ordering and reorder symmetrically
119  ! stored matrix (used when only the nonzero entries in
120  ! the upper triangle of m are being stored)
121  ! 3 reorder symmetrically stored matrix as specified by
122  ! input permutation (used when an ordering has already
123  ! been determined and only the nonzero entries in the
124  ! upper triangle of m are being stored)
125  ! 4 same as 2 but put diagonal entries at start of each row
126  ! 5 same as 3 but put diagonal entries at start of each row
127  !
128  ! flag - integer error flag. values and their meanings are -
129  ! 0 no errors detected
130  ! 9n+k insufficient storage in md
131  ! 10n+1 insufficient storage in odrv
132  ! 11n+1 illegal path specification
133  !
134  !
135  ! conversion from real to double precision
136  !
137  ! change the real declarations in odrv and sro to double precision
138  ! declarations.
139  !
140  !-----------------------------------------------------------------------
141  !
142  implicit none
143 
144  ! -- dummy variables
145  integer(I4B), intent(in) :: n
146  integer(I4B), intent(in) :: nja
147  integer(I4B), intent(in) :: nsp
148  integer(I4B), dimension(n + 1), intent(in) :: ia
149  integer(I4B), dimension(nja), intent(in) :: ja
150  integer(I4B), dimension(n), intent(inout) :: p
151  integer(I4B), dimension(n), intent(inout) :: ip
152  integer(I4B), dimension(nsp), intent(inout) :: isp
153  integer(I4B), intent(inout) :: flag
154 
155  ! -- local
156  integer(I4B) :: v
157  integer(I4B) :: l
158  integer(I4B) :: head
159  integer(I4B) :: mmax
160  integer(I4B) :: next
161  integer(I4B) :: path
162  !
163  ! set path for finding ordering only
164  !
165  path = 1
166  !
167  !
168  ! initialize error flag and validate path specification
169  flag = 0
170  if (path < 1 .or. 5 < path) go to 111
171  !
172  ! find minimum degree ordering
173  mmax = (nsp - n) / 2
174  v = 1
175  l = v + mmax
176  head = l + mmax
177  next = head + n
178  if (mmax < n) go to 110
179  !
180  call ims_md(n, nja, ia, ja, mmax, isp(v), isp(l), isp(head), p, &
181  ip, isp(v), flag)
182  if (flag .ne. 0) go to 100
183  !
184  return
185  !
186  ! ** error -- error detected in md
187  ! flag = 9 * n + vi from routine mdi.
188  !
189 100 return
190  ! ** error -- insufficient storage
191 110 flag = 10 * n + 1
192  return
193  ! ** error -- illegal path specified
194 111 flag = 11 * n + 1
195  return
Here is the call graph for this function:
Here is the caller graph for this function: