Forgot to load lapack in a few examples
[maxima.git] / share / odepack / fortran / md.f
blobfcaaf41d1ca2d949a807cd3ffd3eb87d216d9a37
1 subroutine md
2 * (n, ia,ja, max, v,l, head,last,next, mark, flag)
3 c***********************************************************************
4 c md -- minimum degree algorithm (based on element model)
5 c***********************************************************************
7 c description
9 c md finds a minimum degree ordering of the rows and columns of a
10 c general sparse matrix m stored in (ia,ja,a) format.
11 c when the structure of m is nonsymmetric, the ordering is that
12 c obtained for the symmetric matrix m + m-transpose.
15 c additional parameters
17 c max - declared dimension of the one-dimensional arrays v and l.
18 c max must be at least n+2k, where k is the number of
19 c nonzeroes in the strict upper triangle of m + m-transpose
21 c v - integer one-dimensional work array. dimension = max
23 c l - integer one-dimensional work array. dimension = max
25 c head - integer one-dimensional work array. dimension = n
27 c last - integer one-dimensional array used to return the permutation
28 c of the rows and columns of m corresponding to the minimum
29 c degree ordering. dimension = n
31 c next - integer one-dimensional array used to return the inverse of
32 c the permutation returned in last. dimension = n
34 c mark - integer one-dimensional work array (may be the same as v).
35 c dimension = n
37 c flag - integer error flag. values and their meanings are -
38 c 0 no errors detected
39 c 9n+k insufficient storage in md
42 c definitions of internal parameters
44 c ---------+---------------------------------------------------------
45 c v(s) - value field of list entry
46 c ---------+---------------------------------------------------------
47 c l(s) - link field of list entry (0 =) end of list)
48 c ---------+---------------------------------------------------------
49 c l(vi) - pointer to element list of uneliminated vertex vi
50 c ---------+---------------------------------------------------------
51 c l(ej) - pointer to boundary list of active element ej
52 c ---------+---------------------------------------------------------
53 c head(d) - vj =) vj head of d-list d
54 c - 0 =) no vertex in d-list d
57 c - vi uneliminated vertex
58 c - vi in ek - vi not in ek
59 c ---------+-----------------------------+---------------------------
60 c next(vi) - undefined but nonnegative - vj =) vj next in d-list
61 c - - 0 =) vi tail of d-list
62 c ---------+-----------------------------+---------------------------
63 c last(vi) - (not set until mdp) - -d =) vi head of d-list d
64 c --vk =) compute degree - vj =) vj last in d-list
65 c - ej =) vi prototype of ej - 0 =) vi not in any d-list
66 c - 0 =) do not compute degree -
67 c ---------+-----------------------------+---------------------------
68 c mark(vi) - mark(vk) - nonneg. tag .lt. mark(vk)
71 c - vi eliminated vertex
72 c - ei active element - otherwise
73 c ---------+-----------------------------+---------------------------
74 c next(vi) - -j =) vi was j-th vertex - -j =) vi was j-th vertex
75 c - to be eliminated - to be eliminated
76 c ---------+-----------------------------+---------------------------
77 c last(vi) - m =) size of ei = m - undefined
78 c ---------+-----------------------------+---------------------------
79 c mark(vi) - -m =) overlap count of ei - undefined
80 c - with ek = m -
81 c - otherwise nonnegative tag -
82 c - .lt. mark(vk) -
84 c-----------------------------------------------------------------------
86 integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*),
87 * mark(*), flag, tag, dmin, vk,ek, tail
88 equivalence (vk,ek)
90 c----initialization
91 tag = 0
92 call mdi
93 * (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)
94 if (flag.ne.0) return
96 k = 0
97 dmin = 1
99 c----while k .lt. n do
100 1 if (k.ge.n) go to 4
102 c------search for vertex of minimum degree
103 2 if (head(dmin).gt.0) go to 3
104 dmin = dmin + 1
105 go to 2
107 c------remove vertex vk of minimum degree from degree list
108 3 vk = head(dmin)
109 head(dmin) = next(vk)
110 if (head(dmin).gt.0) last(head(dmin)) = -dmin
112 c------number vertex vk, adjust tag, and tag vk
113 k = k+1
114 next(vk) = -k
115 last(ek) = dmin - 1
116 tag = tag + last(ek)
117 mark(vk) = tag
119 c------form element ek from uneliminated neighbors of vk
120 call mdm
121 * (vk,tail, v,l, last,next, mark)
123 c------purge inactive elements and do mass elimination
124 call mdp
125 * (k,ek,tail, v,l, head,last,next, mark)
127 c------update degrees of uneliminated vertices in ek
128 call mdu
129 * (ek,dmin, v,l, head,last,next, mark)
131 go to 1
133 c----generate inverse permutation from permutation
134 4 do 5 k=1,n
135 next(k) = -next(k)
136 5 last(next(k)) = k
138 return