MODFLOW 6  version 6.5.0.dev2
MODFLOW 6 Code Documentation
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) :: vk
304  integer(I4B) :: ek
305  integer(I4B) :: tail
306  integer(I4B) :: k
307 
308  equivalence(vk, ek)
309  !
310  ! initialization
311  tag = 0
312  call ims_mdi(n, nja, ia, ja, mmax, v, l, head, last, next, &
313  mark, tag, flag)
314  if (flag .ne. 0) return
315  !
316  k = 0
317  dmin = 1
318  !
319  ! while k .lt. n do
320 1 if (k >= n) go to 4
321  !
322  ! search for vertex of minimum degree
323 2 if (head(dmin) > 0) go to 3
324  dmin = dmin + 1
325  go to 2
326  !
327  ! remove vertex vk of minimum degree from degree list
328 3 vk = head(dmin)
329  head(dmin) = next(vk)
330  if (head(dmin) > 0) last(head(dmin)) = -dmin
331  !
332  ! number vertex vk, adjust tag, and tag vk
333  k = k + 1
334  next(vk) = -k
335  last(ek) = dmin - 1
336  tag = tag + last(ek)
337  mark(vk) = tag
338  !
339  ! form element ek from uneliminated neighbors of vk
340  call ims_mdm(n, mmax, vk, tail, v, l, last, next, mark)
341  !
342  ! purge inactive elements and do mass elimination
343  call ims_mdp(n, mmax, k, ek, tail, v, l, head, last, next, mark)
344  !
345  ! update degrees of uneliminated vertices in ek
346  call ims_mdu(n, mmax, ek, dmin, v, l, head, last, next, mark)
347  !
348  go to 1
349  !
350  ! generate inverse permutation from permutation
351 4 do k = 1, n
352  next(k) = -next(k)
353  last(next(k)) = k
354  end do
355  !
356  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 359 of file ImsReordering.f90.

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

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

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

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