1 /* lsodes2.c
2
3 Copyright (c) 1993-2017 Free Software Foundation, Inc.
4
5 This file is part of GNU MCSim.
6
7 GNU MCSim is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation; either version 3
10 of the License, or (at your option) any later version.
11
12 GNU MCSim is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
19
20 lsodes.c was translated from lsodes.f by the utility f2c.
21 To make lsodes.c a stand-alone C routine, the following modifications were
22 made:
23 1. the options -lF77 -lI77 were removed from the link command line
24 2. a function d_sign was written and added at the beginning of
25 the function body
26 3. lsodes was split into two files
27
28 This is the second file for the two parts
29 */
30
31 #include "lsodes.h"
32
33
34 /* ----------------------------------------------------------------------------
35 nsfc_
36
37 symbolic ldu-factorization of nonsymmetric sparse matrix
38 (compressed pointer storage)
39
40 input variables.. n, r, ic, ia, ja, jlmax, jumax.
41 output variables.. il, jl, ijl, iu, ju, iju, flag.
42
43 parameters used internally..
44 nia - q - suppose m* is the result of reordering m. if
45 - processing of the ith row of m* (hence the ith
46 - row of u) is being done, q(j) is initially
47 - nonzero if m*(i,j) is nonzero (j.ge.i). since
48 - values need not be stored, each entry points to the
49 - next nonzero and q(n+1) points to the first. n+1
50 - indicates the end of the list. for example, if n=9
51 - and the 5th row of m* is
52 - 0 x x 0 x 0 0 x 0
53 - then q will initially be
54 - a a a a 8 a a 10 5 (a - arbitrary).
55 - as the algorithm proceeds, other elements of q
56 - are inserted in the list because of fillin.
57 - q is used in an analogous manner to compute the
58 - ith column of l.
59 - size = n+1.
60 nia - ira, - vectors used to find the columns of m. at the kth
61 nia - jra, step of the factorization, irac(k) points to the
62 nia - irac head of a linked list in jra of row indices i
63 - such that i .ge. k and m(i,k) is nonzero. zero
64 - indicates the end of the list. ira(i) (i.ge.k)
65 - points to the smallest j such that j .ge. k and
66 - m(i,j) is nonzero.
67 - size of each = n.
68 nia - irl, - vectors used to find the rows of l. at the kth step
69 nia - jrl of the factorization, jrl(k) points to the head
70 - of a linked list in jrl of column indices j
71 - such j .lt. k and l(k,j) is nonzero. zero
72 - indicates the end of the list. irl(j) (j.lt.k)
73 - points to the smallest i such that i .ge. k and
74 - l(i,j) is nonzero.
75 - size of each = n.
76 nia - iru, - vectors used in a manner analogous to irl and jrl
77 nia - jru to find the columns of u.
78 - size of each = n.
79
80 internal variables..
81 jlptr - points to the last position used in jl.
82 juptr - points to the last position used in ju.
83 jmin,jmax - are the indices in a or u of the first and last
84 elements to be examined in a given row.
85 for example, jmin=ia(k), jmax=ia(k+1)-1.
86 */
87
nsfc_(long * n,long * r,long * ic,long * ia,long * ja,long * jlmax,long * il,long * jl,long * ijl,long * jumax,long * iu,long * ju,long * iju,long * q,long * ira,long * jra,long * irac,long * irl,long * jrl,long * iru,long * jru,long * flag_)88 int nsfc_(long *n, long *r, long *ic, long *ia,
89 long *ja, long *jlmax, long *il, long *jl,
90 long *ijl, long *jumax, long *iu, long *ju,
91 long *iju, long *q, long *ira, long *jra,
92 long *irac, long *irl, long *jrl, long *iru,
93 long *jru, long *flag_)
94 {
95 /* System generated locals */
96 long i__2, i__3;
97
98 /* Local variables */
99 long cend, irai, rend, jmin, jmax, long_, irll, jtmp, irul, i, j, k, m,
100 jaiak, jlmin, lasti, jumin, i1, jlptr, juptr, qm, vj, jairai,
101 lastid, np1, iak, luk;
102 long rk = 0;
103
104 /* Parameter adjustments */
105 --jru;
106 --iru;
107 --jrl;
108 --irl;
109 --irac;
110 --jra;
111 --ira;
112 --q;
113 --iju;
114 --ju;
115 --iu;
116 --ijl;
117 --jl;
118 --il;
119 --ja;
120 --ia;
121 --ic;
122 --r;
123
124 /* Function Body */
125 np1 = *n + 1;
126 jlmin = 1;
127 jlptr = 0;
128 il[1] = 1;
129 jumin = 1;
130 juptr = 0;
131 iu[1] = 1;
132
133 for (k = 1; k <= *n; ++k) {
134 irac[k] = 0;
135 jra[k] = 0;
136 jrl[k] = 0;
137 jru[k] = 0;
138 }
139
140 /* initialize column pointers for a */
141 for (k = 1; k <= *n; ++k) {
142 rk = r[k];
143 iak = ia[rk];
144 if (iak >= ia[rk + 1]) goto L101;
145
146 jaiak = ic[ja[iak]];
147 if (jaiak > k) goto L105;
148
149 jra[k] = irac[jaiak];
150 irac[jaiak] = k;
151
152 ira[k] = iak;
153 }
154
155 /* for each column of l and row of u */
156 for (k = 1; k <= *n; ++k) {
157
158 /* initialize q for computing kth column of l */
159 q[np1] = np1;
160 luk = -1;
161
162 /* by filling in kth column of a */
163 vj = irac[k];
164 if (vj == 0) goto L5;
165
166 L3:
167 qm = np1;
168
169 L4:
170 m = qm;
171 qm = q[m];
172 if (qm < vj) goto L4;
173
174 if (qm == vj) goto L102;
175
176 ++luk;
177 q[m] = vj;
178 q[vj] = qm;
179 vj = jra[vj];
180 if (vj != 0) goto L3;
181
182 /* link through jru */
183 L5:
184 lastid = 0;
185 lasti = 0;
186 ijl[k] = jlptr;
187 i = k;
188
189 L6:
190 i = jru[i];
191 if (i == 0) goto L10;
192
193 qm = np1;
194 jmin = irl[i];
195 jmax = ijl[i] + il[i + 1] - il[i] - 1;
196 long_ = jmax - jmin;
197 if (long_ < 0) goto L6;
198
199 jtmp = jl[jmin];
200 if (jtmp != k) ++long_;
201
202 if (jtmp == k) r[i] = -r[i];
203
204 if (lastid < long_) {
205 lasti = i;
206 lastid = long_;
207 }
208
209 /* and merge the corresponding columns into the kth column */
210 for (j = jmin; j <= jmax; ++j) {
211 vj = jl[j];
212
213 do {
214 m = qm;
215 qm = q[m];
216 }
217 while (qm < vj);
218
219 if (qm != vj) {
220 ++luk;
221 q[m] = vj;
222 q[vj] = qm;
223 qm = vj;
224 }
225 }
226
227 goto L6;
228
229 L10:
230 /* lasti is the longest column merged into the kth.
231 see if it equals the entire kth column */
232 qm = q[np1];
233 if (qm != k) goto L105;
234
235 if (luk == 0) goto L17;
236
237 if (lastid != luk) goto L11;
238
239 /* if so, jl can be compressed */
240 irll = irl[lasti];
241 ijl[k] = irll + 1;
242 if (jl[irll] != k) --ijl[k];
243
244 goto L17;
245
246 L11:
247 /* if not, see if kth column can overlap the previous one */
248 if (jlmin > jlptr) goto L15;
249
250 qm = q[qm];
251 for (j = jlmin; j <= jlptr; ++j) {
252 if ((i__3 = jl[j] - qm) >= 0) {
253 if (i__3 == 0) goto L13;
254 else goto L15;
255 }
256 }
257
258 goto L15;
259
260 L13:
261 ijl[k] = j;
262 for (i = j; i <= jlptr; ++i) {
263 if (jl[i] != qm) goto L15;
264
265 qm = q[qm];
266 if (qm > *n) goto L17;
267 }
268
269 jlptr = j - 1;
270
271 L15:
272 /* * move column indices from q to jl, update vectors **
273 * */
274 jlmin = jlptr + 1;
275 ijl[k] = jlmin;
276 if (luk == 0) {
277 goto L17;
278 }
279 jlptr += luk;
280 if (jlptr > *jlmax) {
281 goto L103;
282 }
283 qm = q[np1];
284 i__2 = jlptr;
285 for (j = jlmin; j <= i__2; ++j) {
286 qm = q[qm];
287 /* L16: */
288 jl[j] = qm;
289 }
290 L17:
291 irl[k] = ijl[k];
292 il[k + 1] = il[k] + luk;
293
294 /* * initialize q for computing kth row of u *
295 * */
296 q[np1] = np1;
297 luk = -1;
298 /* * by filling in kth row of reordered a *
299 * */
300 rk = r[k];
301 jmin = ira[k];
302 jmax = ia[rk + 1] - 1;
303 if (jmin > jmax) {
304 goto L20;
305 }
306 i__2 = jmax;
307 for (j = jmin; j <= i__2; ++j) {
308 vj = ic[ja[j]];
309 qm = np1;
310 L18:
311 m = qm;
312 qm = q[m];
313 if (qm < vj) {
314 goto L18;
315 }
316 if (qm == vj) {
317 goto L102;
318 }
319 ++luk;
320 q[m] = vj;
321 q[vj] = qm;
322 /* L19: */
323 }
324 /* * link through jrl, **
325 * */
326 L20:
327 lastid = 0;
328 lasti = 0;
329 iju[k] = juptr;
330 i = k;
331 i1 = jrl[k];
332 L21:
333 i = i1;
334 if (i == 0) {
335 goto L26;
336 }
337 i1 = jrl[i];
338 qm = np1;
339 jmin = iru[i];
340 jmax = iju[i] + iu[i + 1] - iu[i] - 1;
341 long_ = jmax - jmin;
342 if (long_ < 0) {
343 goto L21;
344 }
345 jtmp = ju[jmin];
346 if (jtmp == k) {
347 goto L22;
348 }
349 /* * update irl and jrl, **
350 * */
351 ++long_;
352 cend = ijl[i] + il[i + 1] - il[i];
353 ++irl[i];
354 if (irl[i] >= cend) {
355 goto L22;
356 }
357 j = jl[irl[i]];
358 jrl[i] = jrl[j];
359 jrl[j] = i;
360 L22:
361 if (lastid >= long_) {
362 goto L23;
363 }
364 lasti = i;
365 lastid = long_;
366 /* * and merge the corresponding rows into the kth row *
367 * */
368 L23:
369 i__2 = jmax;
370 for (j = jmin; j <= i__2; ++j) {
371 vj = ju[j];
372 L24:
373 m = qm;
374 qm = q[m];
375 if (qm < vj) {
376 goto L24;
377 }
378 if (qm == vj) {
379 goto L25;
380 }
381 ++luk;
382 q[m] = vj;
383 q[vj] = qm;
384 qm = vj;
385 L25:
386 ;
387 }
388 goto L21;
389 /* * update jrl(k) and irl(k) *
390 * */
391 L26:
392 if (il[k + 1] <= il[k]) {
393 goto L27;
394 }
395 j = jl[irl[k]];
396 jrl[k] = jrl[j];
397 jrl[j] = k;
398 /* * lasti is the longest row merged into the kth *
399 * */
400 /* * see if it equals the entire kth row **
401 * */
402 L27:
403 qm = q[np1];
404 if (qm != k) {
405 goto L105;
406 }
407 if (luk == 0) {
408 goto L34;
409 }
410 if (lastid != luk) {
411 goto L28;
412 }
413 /* * if so, ju can be compressed **
414 * */
415 irul = iru[lasti];
416 iju[k] = irul + 1;
417 if (ju[irul] != k) {
418 --iju[k];
419 }
420 goto L34;
421
422 L28:
423 /* if not, see if kth row can overlap the previous one */
424 if (jumin > juptr) goto L32;
425
426 qm = q[qm];
427 for (j = jumin; j <= juptr; ++j) {
428 if ((i__3 = ju[j] - qm) < 0) goto L29;
429 else
430 if (i__3 == 0) goto L30;
431 else goto L32;
432
433 L29: ;
434 }
435 goto L32;
436
437 L30:
438 iju[k] = j;
439 for (i = j; i <= juptr; ++i) {
440 if (ju[i] != qm) goto L32;
441
442 qm = q[qm];
443 if (qm > *n) goto L34;
444 }
445 juptr = j - 1;
446
447 L32:
448 /* move row indices from q to ju, update vectors */
449 jumin = juptr + 1;
450 iju[k] = jumin;
451 if (luk == 0) goto L34;
452
453 juptr += luk;
454 if (juptr > *jumax) goto L106;
455
456 qm = q[np1];
457 for (j = jumin; j <= juptr; ++j) {
458 qm = q[qm];
459 ju[j] = qm;
460 }
461
462 L34:
463 iru[k] = iju[k];
464 iu[k + 1] = iu[k] + luk;
465
466 /* update iru, jru */
467 i = k;
468
469 L35:
470 i1 = jru[i];
471 if (r[i] < 0) goto L36;
472
473 rend = iju[i] + iu[i + 1] - iu[i];
474 if (iru[i] >= rend) goto L37;
475
476 j = ju[iru[i]];
477 jru[i] = jru[j];
478 jru[j] = i;
479 goto L37;
480
481 L36:
482 r[i] = -r[i];
483
484 L37:
485 i = i1;
486 if (i == 0) goto L38;
487
488 ++iru[i];
489 goto L35;
490
491 L38:
492 /* update ira, jra, irac */
493 i = irac[k];
494 if (i == 0) goto L41;
495
496 L39:
497 i1 = jra[i];
498 ++ira[i];
499 if (ira[i] >= ia[r[i] + 1]) goto L40;
500
501 irai = ira[i];
502 jairai = ic[ja[irai]];
503 if (jairai > i) goto L40;
504
505 jra[i] = irac[jairai];
506 irac[jairai] = i;
507
508 L40:
509 i = i1;
510 if (i != 0) goto L39;
511
512 L41: ;
513 }
514
515 ijl[*n] = jlptr;
516 iju[*n] = juptr;
517 *flag_ = 0;
518 return 0;
519
520 L101:
521 /* error.. null row in a */
522 *flag_ = *n + rk;
523 return 0;
524
525 L102:
526 /* error.. duplicate entry in a */
527 *flag_ = (*n << 1) + rk;
528 return 0;
529
530 L103:
531 /* error.. insufficient storage for jl */
532 *flag_ = *n * 3 + k;
533 return 0;
534
535 L105:
536 /* error.. null pivot */
537 *flag_ = *n * 5 + k;
538 return 0;
539
540 L106:
541 /* error.. insufficient storage for ju */
542 *flag_ = *n * 6 + k;
543 return 0;
544 } /* nsfc_ */
545
nroc_(long * n,long * ic,long * ia,long * ja,double * a,long * jar,double * ar,long * p,long * flag_)546 int nroc_(long *n, long *ic, long *ia, long *ja,
547 double *a, long *jar, double *ar, long *p,
548 long *flag_)
549 {
550 /* System generated locals */
551 long i__1, i__2;
552
553 /* Local variables */
554 long jmin, jmax, newj, i, j, k;
555
556 /* -----------------------------------------------------------------------------
557 */
558 /* 0 0 0 0 d */
559 /* would be represented as */
560 /* - 1 2 3 4 5 6 */
561 /* ----+------------ */
562 /* iu - 1 4 6 7 8 8 */
563 /* ju - 3 4 5 4 */
564 /* iju - 1 2 4 3 . */
565 /* the diagonal entries of l and u are assumed to be equal to one and
566 */
567 /* are not stored. the array d contains the reciprocals of the */
568 /* diagonal entries of the matrix d. */
569
570 /* iii. additional storage savings */
571 /* in nsfc, r and ic can be the same array in the calling */
572 /* sequence if no reordering of the coefficient matrix has been done.
573 */
574 /* in nnfc, r, c, and ic can all be the same array if no */
575 /* reordering has been done. if only the rows have been reordered, */
576 /* then c and ic can be the same array. if the row and column */
577 /* orderings are the same, then r and c can be the same array. z and
578 */
579 /* row can be the same array. */
580 /* in nnsc or nntc, r and c can be the same array if no */
581 /* reordering has been done or if the row and column orderings are the
582 */
583 /* same. z and b can be the same array. however, then b will be */
584 /* destroyed. */
585
586 /* iv. parameters */
587 /* following is a list of parameters to the programs. names are
588 */
589 /* uniform among the various subroutines. class abbreviations are --
590 */
591 /* n - long variable */
592 /* f - real variable */
593 /* v - supplies a value to a subroutine */
594 /* r - returns a result from a subroutine */
595 /* i - used internally by a subroutine */
596 /* a - array */
597
598 /* class - parameter */
599 /* ------+---------- */
600 /* fva - a - nonzero entries of the coefficient matrix m, stored */
601 /* - by rows. */
602 /* - size = number of nonzero entries in m. */
603 /* fva - b - right-hand side b. */
604 /* - size = n. */
605 /* nva - c - ordering of the columns of m. */
606 /* - size = n. */
607 /* fvra - d - reciprocals of the diagonal entries of the matrix d. */
608
609 /* - size = n. */
610 /* nr - flag - error flag. values and their meanings are -- */
611 /* - 0 no errors detected */
612 /* - n+k null row in a -- row = k */
613 /* - 2n+k duplicate entry in a -- row = k */
614 /* - 3n+k insufficient storage for jl -- row = k */
615 /* - 4n+1 insufficient storage for l */
616 /* - 5n+k null pivot -- row = k */
617 /* - 6n+k insufficient storage for ju -- row = k */
618 /* - 7n+1 insufficient storage for u */
619 /* - 8n+k zero pivot -- row = k */
620 /* nva - ia - pointers to delimit the rows of a. */
621 /* - size = n+1. */
622 /* nvra - ijl - pointers to the first element in each column in jl, */
623 /* - used to compress storage in jl. */
624 /* - size = n. */
625 /* nvra - iju - pointers to the first element in each row in ju, used
626 */
627 /* - to compress storage in ju. */
628 /* - size = n. */
629 /* nvra - il - pointers to delimit the columns of l. */
630 /* - size = n+1. */
631 /* nvra - iu - pointers to delimit the rows of u. */
632 /* - size = n+1. */
633 /* nva - ja - column numbers corresponding to the elements of a. */
634 /* - size = size of a. */
635 /* nvra - jl - row numbers corresponding to the elements of l. */
636 /* - size = jlmax. */
637 /* nv - jlmax - declared dimension of jl. jlmax must be larger than */
638
639 /* - the number of nonzeros in the strict lower triangle
640 */
641 /* - of m plus fillin minus compression. */
642 /* nvra - ju - column numbers corresponding to the elements of u. */
643 /* - size = jumax. */
644 /* nv - jumax - declared dimension of ju. jumax must be larger than */
645
646 /* - the number of nonzeros in the strict upper triangle
647 */
648 /* - of m plus fillin minus compression. */
649 /* fvra - l - nonzero entries in the strict lower triangular portion
650 */
651 /* - of the matrix l, stored by columns. */
652 /* - size = lmax. */
653 /* nv - lmax - declared dimension of l. lmax must be larger than */
654 /* - the number of nonzeros in the strict lower triangle
655 */
656 /* - of m plus fillin (il(n+1)-1 after nsfc). */
657 /* nv - n - number of variables/equations. */
658 /* nva - r - ordering of the rows of m. */
659 /* - size = n. */
660 /* fvra - u - nonzero entries in the strict upper triangular portion
661 */
662 /* - of the matrix u, stored by rows. */
663 /* - size = umax. */
664 /* nv - umax - declared dimension of u. umax must be larger than */
665 /* - the number of nonzeros in the strict upper triangle
666 */
667 /* - of m plus fillin (iu(n+1)-1 after nsfc). */
668 /* fra - z - solution x. */
669 /* - size = n. */
670
671 /* ----------------------------------------------------------------
672 */
673
674 /* subroutine nroc */
675 /* reorders rows of a, leaving row order unchanged */
676
677
678 /* input parameters.. n, ic, ia, ja, a */
679 /* output parameters.. ja, a, flag */
680
681 /* parameters used internally.. */
682 /* nia - p - at the kth step, p is a linked list of the reordered */
683
684 /* - column indices of the kth row of a. p(n+1) points */
685
686 /* - to the first entry in the list. */
687 /* - size = n+1. */
688 /* nia - jar - at the kth step,jar contains the elements of the */
689 /* - reordered column indices of a. */
690 /* - size = n. */
691 /* fia - ar - at the kth step, ar contains the elements of the */
692 /* - reordered row of a. */
693 /* - size = n. */
694
695
696 /* * for each nonempty row * */
697 /* Parameter adjustments */
698 --p;
699 --ar;
700 --jar;
701 --a;
702 --ja;
703 --ia;
704 --ic;
705
706 /* Function Body */
707 i__1 = *n;
708 for (k = 1; k <= i__1; ++k) {
709 jmin = ia[k];
710 jmax = ia[k + 1] - 1;
711 if (jmin > jmax) {
712 goto L5;
713 }
714 p[*n + 1] = *n + 1;
715 /* * insert each element in the list * */
716 i__2 = jmax;
717 for (j = jmin; j <= i__2; ++j) {
718 newj = ic[ja[j]];
719 i = *n + 1;
720 L1:
721 if (p[i] >= newj) {
722 goto L2;
723 }
724 i = p[i];
725 goto L1;
726 L2:
727 if (p[i] == newj) {
728 goto L102;
729 }
730 p[newj] = p[i];
731 p[i] = newj;
732 jar[newj] = ja[j];
733 ar[newj] = a[j];
734 /* L3: */
735 }
736 /* * replace old row in ja and a ** */
737 i = *n + 1;
738 i__2 = jmax;
739 for (j = jmin; j <= i__2; ++j) {
740 i = p[i];
741 ja[j] = jar[i];
742 /* L4: */
743 a[j] = ar[i];
744 }
745 L5:
746 ;
747 }
748 *flag_ = 0;
749 return 0;
750
751 /* error.. duplicate entry in a */
752 L102:
753 *flag_ = *n + k;
754 return 0;
755 } /* nroc_ */
756
757 /* -----------------------------------------------------------------------------
758 odrv_ 5/2/83
759 */
760 /* the permutation returned in p. dimension = n */
761
762 /* nsp - declared dimension of the one-dimensional array isp. nsp */
763 /* must be at least 3n+4k, where k is the number of nonzeroes
764 */
765 /* in the strict upper triangle of m */
766
767 /* isp - long one-dimensional array used for working storage. */
768 /* dimension = nsp */
769
770 /* path - long path specification. values and their meanings are -
771 */
772 /* 1 find minimum degree ordering only */
773 /* 2 find minimum degree ordering and reorder symmetrically
774 */
775 /* stored matrix (used when only the nonzero entries in
776 */
777 /* the upper triangle of m are being stored) */
778 /* 3 reorder symmetrically stored matrix as specified by */
779 /* input permutation (used when an ordering has already
780 */
781 /* been determined and only the nonzero entries in the */
782
783 /* upper triangle of m are being stored) */
784 /* 4 same as 2 but put diagonal entries at start of each row
785 */
786 /* 5 same as 3 but put diagonal entries at start of each row
787 */
788
789 /* flag - long error flag. values and their meanings are - */
790 /* 0 no errors detected */
791 /* 9n+k insufficient storage in md */
792 /* 10n+1 insufficient storage in odrv */
793 /* 11n+1 illegal path specification */
794
795
796 /* conversion from real to double precision */
797
798 /* change the real declarations in odrv and sro to double precision */
799 /* declarations. */
800
801 /* -------------------------------------------------------------------
802 */
odrv_(long * n,long * ia,long * ja,double * a,long * p,long * ip,long * nsp,long * isp,long path,long * flag_)803 int odrv_(long *n, long *ia, long *ja, double *a,
804 long *p, long *ip, long *nsp, long *isp,
805 long path, long *flag_)
806 {
807 long head, l;
808 long dflag;
809 long q, v;
810 extern /* Subroutine */ int md_(long *, long *, long *, long *
811 , long *, long *, long *, long *, long *, long *
812 , long *);
813 long max_, tmp;
814 extern /* Subroutine */ int sro_(long *, long *, long *, long
815 *, double *, long *, long *, long *);
816
817
818 /* ----initialize error flag and validate path specification */
819 /* Parameter adjustments */
820 --isp;
821 --ip;
822 --p;
823 --a;
824 --ja;
825 --ia;
826
827 /* Function Body */
828 *flag_ = 0;
829 if (path < 1 || 5 < path) {
830 goto L111;
831 }
832
833 /* ----allocate storage and find minimum degree ordering */
834 if ((path - 1) * (path - 2) * (path - 4) != 0) {
835 goto L1;
836 }
837 max_ = (*nsp - *n) / 2;
838 v = 1;
839 l = v + max_;
840 head = l + max_;
841 if (max_ < *n) goto L110;
842
843 md_(n, &ia[1], &ja[1], &max_, &isp[v], &isp[l], &isp[head], &p[1], &ip[1],
844 &isp[v], flag_);
845 if (*flag_ != 0) goto L100;
846
847 /* ----allocate storage and symmetrically reorder matrix */
848 L1:
849 if ((path - 2) * (path - 3) * (path - 4) * (path - 5) != 0) {
850 goto L2;
851 }
852 tmp = *nsp + 1 - *n;
853 q = tmp - (ia[*n + 1] - 1);
854 if (q < 1) {
855 goto L110;
856 }
857
858 dflag = path == 4 || path == 5;
859 sro_(n, &ip[1], &ia[1], &ja[1], &a[1], &isp[tmp], &isp[q], &dflag);
860
861 L2:
862 return 0;
863
864 /* error -- error detected in md */
865 L100:
866 return 0;
867 /* error -- insufficient storage */
868 L110:
869 *flag_ = *n * 10 + 1;
870 return 0;
871 /* error -- illegal path specified */
872 L111:
873 *flag_ = *n * 11 + 1;
874 return 0;
875 } /* odrv_ */
876
sro_(long * n,long * ip,long * ia,long * ja,double * a,long * q,long * r,long * dflag)877 int sro_(long *n, long *ip, long *ia, long *ja,
878 double *a, long *q, long *r, long *dflag)
879 {
880 /* System generated locals */
881 long i__1, i__2;
882
883 /* Local variables */
884 long jmin, jmax, i, j, k, ilast;
885 double ak;
886 long jdummy, jak;
887
888 /*
889 */
890 /* sro -- symmetric reordering of sparse symmetric matrix */
891 /*
892 */
893
894 /* description */
895
896 /* the nonzero entries of the matrix m are assumed to be stored */
897 /* symmetrically in (ia,ja,a) format (i.e., not both m(i,j) and m(j,i)
898 */
899 /* are stored if i ne j). */
900
901 /* sro does not rearrange the order of the rows, but does move */
902 /* nonzeroes from one row to another to ensure that if m(i,j) will be
903 */
904 /* in the upper triangle of m with respect to the new ordering, then */
905
906 /* m(i,j) is stored in row i (and thus m(j,i) is not stored), whereas
907 */
908 /* if m(i,j) will be in the strict lower triangle of m, then m(j,i) is
909 */
910 /* stored in row j (and thus m(i,j) is not stored). */
911
912
913 /* additional parameters */
914
915 /* q - long one-dimensional work array. dimension = n */
916
917 /* r - long one-dimensional work array. dimension = number of
918 */
919 /* nonzero entries in the upper triangle of m */
920
921 /* dflag - logical variable. if dflag = .true., then store nonzero */
922 /* diagonal elements at the beginning of the row */
923
924 /* -------------------------------------------------------------------
925 */
926
927
928 /* --phase 1 -- find row in which to store each nonzero */
929 /* ----initialize count of nonzeroes to be stored in each row */
930 /* Parameter adjustments */
931 --r;
932 --q;
933 --a;
934 --ja;
935 --ia;
936 --ip;
937
938 /* Function Body */
939 i__1 = *n;
940 for (i = 1; i <= i__1; ++i) {
941 /* L1: */
942 q[i] = 0;
943 }
944
945 /* ----for each nonzero element a(j) */
946 i__1 = *n;
947 for (i = 1; i <= i__1; ++i) {
948 jmin = ia[i];
949 jmax = ia[i + 1] - 1;
950 if (jmin > jmax) {
951 goto L3;
952 }
953 i__2 = jmax;
954 for (j = jmin; j <= i__2; ++j) {
955
956 /* --------find row (=r(j)) and column (=ja(j)) in which to store
957 a(j) ... */
958 k = ja[j];
959 if (ip[k] < ip[i]) {
960 ja[j] = i;
961 }
962 if (ip[k] >= ip[i]) {
963 k = i;
964 }
965 r[j] = k;
966
967 /* --------... and increment count of nonzeroes (=q(r(j)) in that
968 row */
969 /* L2: */
970 ++q[k];
971 }
972 L3:
973 ;
974 }
975
976
977 /* --phase 2 -- find new ia and permutation to apply to (ja,a) */
978 /* ----determine pointers to delimit rows in permuted (ja,a) */
979 i__1 = *n;
980 for (i = 1; i <= i__1; ++i) {
981 ia[i + 1] = ia[i] + q[i];
982 /* L4: */
983 q[i] = ia[i + 1];
984 }
985
986 /* ----determine where each (ja(j),a(j)) is stored in permuted (ja,a) */
987 /* ----for each nonzero element (in reverse order) */
988 ilast = 0;
989 jmin = ia[1];
990 jmax = ia[*n + 1] - 1;
991 j = jmax;
992 i__1 = jmax;
993 for (jdummy = jmin; jdummy <= i__1; ++jdummy) {
994 i = r[j];
995 if (! (*dflag) || ja[j] != i || i == ilast) {
996 goto L5;
997 }
998
999 /* ------if dflag, then put diagonal nonzero at beginning of row */
1000 r[j] = ia[i];
1001 ilast = i;
1002 goto L6;
1003
1004 /* ------put (off-diagonal) nonzero in last unused location in row */
1005 L5:
1006 --q[i];
1007 r[j] = q[i];
1008
1009 L6:
1010 --j;
1011 }
1012
1013
1014 /* --phase 3 -- permute (ja,a) to upper triangular form (wrt new ordering)
1015 */
1016 i__1 = jmax;
1017 for (j = jmin; j <= i__1; ++j) {
1018 L7:
1019 if (r[j] == j) {
1020 goto L8;
1021 }
1022 k = r[j];
1023 r[j] = r[k];
1024 r[k] = k;
1025 jak = ja[k];
1026 ja[k] = ja[j];
1027 ja[j] = jak;
1028 ak = a[k];
1029 a[k] = a[j];
1030 a[j] = ak;
1031 goto L7;
1032 L8:
1033 ;
1034 }
1035
1036 return 0;
1037 } /* sro_ */
1038
md_(long * n,long * ia,long * ja,long * max_,long * v,long * l,long * head,long * last,long * next,long * mark,long * flag_)1039 int md_(long *n, long *ia, long *ja, long *max_,
1040 long *v, long *l, long *head, long *last,
1041 long *next, long *mark, long *flag_)
1042 {
1043 /* System generated locals */
1044 long i__1;
1045 static long equiv_0[1];
1046
1047 /* Local variables */
1048 long dmin_, tail, k;
1049 #define ek (equiv_0)
1050 #define vk (equiv_0)
1051 long tag;
1052
1053 /*
1054 */
1055 /* md -- minimum degree algorithm (based on element model) */
1056 /*
1057 */
1058
1059 /* description */
1060
1061 /* md finds a minimum degree ordering of the rows and columns of a */
1062 /* general sparse matrix m stored in (ia,ja,a) format. */
1063 /* when the structure of m is nonsymmetric, the ordering is that */
1064 /* obtained for the symmetric matrix m + m-transpose. */
1065
1066
1067 /* additional parameters */
1068
1069 /* max - declared dimension of the one-dimensional arrays v and l. */
1070 /* max must be at least n+2k, where k is the number of */
1071 /* nonzeroes in the strict upper triangle of m + m-transpose */
1072
1073 /* v - long one-dimensional work array. dimension = max */
1074
1075 /* l - long one-dimensional work array. dimension = max */
1076
1077 /* head - long one-dimensional work array. dimension = n */
1078
1079 /* last - long one-dimensional array used to return the permutation
1080 */
1081 /* of the rows and columns of m corresponding to the minimum */
1082 /* degree ordering. dimension = n */
1083
1084 /* next - long one-dimensional array used to return the inverse of
1085 */
1086 /* the permutation returned in last. dimension = n */
1087
1088 /* mark - long one-dimensional work array (may be the same as v). */
1089
1090 /* dimension = n */
1091
1092 /* flag - long error flag. values and their meanings are - */
1093 /* 0 no errors detected */
1094 /* 9n+k insufficient storage in md */
1095
1096
1097 /* definitions of internal parameters */
1098
1099 /* ---------+---------------------------------------------------------
1100 */
1101 /* v(s) - value field of list entry */
1102 /* ---------+---------------------------------------------------------
1103 */
1104 /* l(s) - link field of list entry (0 =) end of list) */
1105 /* ---------+---------------------------------------------------------
1106 */
1107 /* l(vi) - pointer to element list of uneliminated vertex vi */
1108 /* ---------+---------------------------------------------------------
1109 */
1110 /* l(ej) - pointer to boundary list of active element ej */
1111 /* ---------+---------------------------------------------------------
1112 */
1113 /* head(d) - vj =) vj head of d-list d */
1114 /* - 0 =) no vertex in d-list d */
1115
1116
1117 /* - vi uneliminated vertex */
1118 /* - vi in ek - vi not in ek */
1119 /* ---------+-----------------------------+---------------------------
1120 */
1121 /* next(vi) - undefined but nonnegative - vj =) vj next in d-list */
1122 /* - - 0 =) vi tail of d-list */
1123 /* ---------+-----------------------------+---------------------------
1124 */
1125 /* last(vi) - (not set until mdp) - -d =) vi head of d-list d
1126 */
1127 /* --vk =) compute degree - vj =) vj last in d-list */
1128 /* - ej =) vi prototype of ej - 0 =) vi not in any d-list
1129 */
1130 /* - 0 =) do not compute degree - */
1131 /* ---------+-----------------------------+---------------------------
1132 */
1133 /* mark(vi) - mark(vk) - nonneg. tag .lt. mark(vk)
1134 */
1135
1136
1137 /* - vi eliminated vertex */
1138 /* - ei active element - otherwise */
1139 /* ---------+-----------------------------+---------------------------
1140 */
1141 /* next(vi) - -j =) vi was j-th vertex - -j =) vi was j-th vertex */
1142
1143 /* - to be eliminated - to be eliminated */
1144 /* ---------+-----------------------------+---------------------------
1145 */
1146 /* last(vi) - m =) size of ei = m - undefined */
1147 /* ---------+-----------------------------+---------------------------
1148 */
1149 /* mark(vi) - -m =) overlap count of ei - undefined */
1150 /* - with ek = m - */
1151 /* - otherwise nonnegative tag - */
1152 /* - .lt. mark(vk) - */
1153
1154 /* -------------------------------------------------------------------
1155 */
1156
1157
1158 /* ----initialization */
1159 /* Parameter adjustments */
1160 --mark;
1161 --next;
1162 --last;
1163 --head;
1164 --l;
1165 --v;
1166 --ja;
1167 --ia;
1168
1169 /* Function Body */
1170 tag = 0;
1171 mdi_(n, &ia[1], &ja[1], max_, &v[1], &l[1], &head[1], &last[1], &next[1],
1172 &mark[1], &tag, flag_);
1173 if (*flag_ != 0) {
1174 return 0;
1175 }
1176
1177 k = 0;
1178 dmin_ = 1;
1179
1180 /* while k .lt. n do */
1181 L1:
1182 if (k >= *n) {
1183 goto L4;
1184 }
1185
1186 /* search for vertex of minimum degree */
1187 L2:
1188 if (head[dmin_] > 0) {
1189 goto L3;
1190 }
1191 ++dmin_;
1192 goto L2;
1193
1194 /* remove vertex vk of minimum degree from degree list */
1195 L3:
1196 *vk = head[dmin_];
1197 head[dmin_] = next[*vk];
1198 if (head[dmin_] > 0) {
1199 last[head[dmin_]] = -dmin_;
1200 }
1201
1202 /* ------number vertex vk, adjust tag, and tag vk */
1203 ++k;
1204 next[*vk] = -k;
1205 last[*ek] = dmin_ - 1;
1206 tag += last[*ek];
1207 mark[*vk] = tag;
1208
1209 /* ------form element ek from uneliminated neighbors of vk */
1210 mdm_(vk, &tail, &v[1], &l[1], &last[1], &next[1], &mark[1]);
1211
1212 /* ------purge inactive elements and do mass elimination */
1213 mdp_(&k, ek, &tail, &v[1], &l[1], &head[1], &last[1], &next[1], &mark[1]);
1214
1215
1216 /* ------update degrees of uneliminated vertices in ek */
1217 mdu_(ek, &dmin_, &v[1], &l[1], &head[1], &last[1], &next[1], &mark[1]);
1218
1219 goto L1;
1220
1221 /* ----generate inverse permutation from permutation */
1222 L4:
1223 i__1 = *n;
1224 for (k = 1; k <= i__1; ++k) {
1225 next[k] = -next[k];
1226 /* L5: */
1227 last[next[k]] = k;
1228 }
1229
1230 return 0;
1231 } /* md_ */
1232
1233 #undef vk
1234 #undef ek
1235
1236
mdu_(long * ek,long * dmin_,long * v,long * l,long * head,long * last,long * next,long * mark)1237 int mdu_(long *ek, long *dmin_, long *v, long *l,
1238 long *head, long *last, long *next, long *mark)
1239 {
1240 /* System generated locals */
1241 long i__1, i__2;
1242 static long equiv_0[1];
1243
1244 /* Local variables */
1245 long b, i, s;
1246 #define es (equiv_0)
1247 long vb, vi;
1248 #define vs (equiv_0)
1249 long blpmax, ilpmax, tag, blp, dvi, evi, ilp;
1250
1251 /* lll. optimize */
1252 /*
1253 */
1254 /* mdu -- update degrees of uneliminated vertices in ek */
1255 /*
1256 */
1257
1258 /* ----initialize tag */
1259 /* Parameter adjustments */
1260 --mark;
1261 --next;
1262 --last;
1263 --head;
1264 --l;
1265 --v;
1266
1267 /* Function Body */
1268 tag = mark[*ek] - last[*ek];
1269
1270 /* ----for each vertex vi in ek */
1271 i = *ek;
1272 ilpmax = last[*ek];
1273 if (ilpmax <= 0) {
1274 goto L11;
1275 }
1276 i__1 = ilpmax;
1277 for (ilp = 1; ilp <= i__1; ++ilp) {
1278 i = l[i];
1279 vi = v[i];
1280 if ((i__2 = last[vi]) < 0) {
1281 goto L1;
1282 } else if (i__2 == 0) {
1283 goto L10;
1284 } else goto L8;
1285
1286 /* if vi neither prototype nor duplicate vertex, then merge elements */
1287 /* to compute degree */
1288 L1:
1289 ++tag;
1290 dvi = last[*ek];
1291
1292 /* for each vertex/element vs/es in element list of vi */
1293 s = l[vi];
1294 L2:
1295 s = l[s];
1296 if (s == 0) {
1297 goto L9;
1298 }
1299 *vs = v[s];
1300 if (next[*vs] < 0) {
1301 goto L3;
1302 }
1303
1304 /* if vs is uneliminated vertex, then tag and adjust degree */
1305 mark[*vs] = tag;
1306 ++dvi;
1307 goto L5;
1308
1309 /* ----------if es is active element, then expand */
1310 /* ------------check for outmatched vertex */
1311 L3:
1312 if (mark[*es] < 0) {
1313 goto L6;
1314 }
1315
1316 /* ------------for each vertex vb in es */
1317 b = *es;
1318 blpmax = last[*es];
1319 i__2 = blpmax;
1320 for (blp = 1; blp <= i__2; ++blp) {
1321 b = l[b];
1322 vb = v[b];
1323
1324 /* --------------if vb is untagged, then tag and adjust degree */
1325 if (mark[vb] >= tag) {
1326 goto L4;
1327 }
1328 mark[vb] = tag;
1329 ++dvi;
1330 L4:
1331 ;
1332 }
1333
1334 L5:
1335 goto L2;
1336
1337 /* ------else if vi is outmatched vertex, then adjust overlaps but do
1338 not */
1339 /* ------compute degree */
1340 L6:
1341 last[vi] = 0;
1342 --mark[*es];
1343 L7:
1344 s = l[s];
1345 if (s == 0) {
1346 goto L10;
1347 }
1348 *es = v[s];
1349 if (mark[*es] < 0) {
1350 --mark[*es];
1351 }
1352 goto L7;
1353
1354 /* ------else if vi is prototype vertex, then calculate degree by */
1355 /* ------inclusion/exclusion and reset overlap count */
1356 L8:
1357 evi = last[vi];
1358 dvi = last[*ek] + last[evi] + mark[evi];
1359 mark[evi] = 0;
1360
1361 /* ------insert vi in appropriate degree list */
1362 L9:
1363 next[vi] = head[dvi];
1364 head[dvi] = vi;
1365 last[vi] = -dvi;
1366 if (next[vi] > 0) {
1367 last[next[vi]] = vi;
1368 }
1369 if (dvi < *dmin_) {
1370 *dmin_ = dvi;
1371 }
1372
1373 L10:
1374 ;
1375 }
1376
1377 L11:
1378 return 0;
1379 } /* mdu_ */
1380
1381 #undef vs
1382 #undef es
1383
1384
mdp_(long * k,long * ek,long * tail,long * v,long * l,long * head,long * last,long * next,long * mark)1385 int mdp_(long *k, long *ek, long *tail, long *v,
1386 long *l, long *head, long *last, long *next,
1387 long *mark)
1388 {
1389 /* System generated locals */
1390 long i__1;
1391
1392 /* Local variables */
1393 long i, s, li, es, vi, ls, ilpmax, tag, evi, ilp, lvi;
1394 long free = 0;
1395
1396 /* mdp -- purge inactive elements and do mass elimination */
1397
1398 /* ----initialize tag */
1399 /* Parameter adjustments */
1400 --mark;
1401 --next;
1402 --last;
1403 --head;
1404 --l;
1405 --v;
1406
1407 /* Function Body */
1408 tag = mark[*ek];
1409
1410 /* ----for each vertex vi in ek */
1411 li = *ek;
1412 ilpmax = last[*ek];
1413 if (ilpmax <= 0) {
1414 goto L12;
1415 }
1416 i__1 = ilpmax;
1417 for (ilp = 1; ilp <= i__1; ++ilp) {
1418 i = li;
1419 li = l[i];
1420 vi = v[li];
1421
1422 /* ------remove vi from degree list */
1423 if (last[vi] == 0) {
1424 goto L3;
1425 }
1426 if (last[vi] > 0) {
1427 goto L1;
1428 }
1429 head[-last[vi]] = next[vi];
1430 goto L2;
1431 L1:
1432 next[last[vi]] = next[vi];
1433 L2:
1434 if (next[vi] > 0) {
1435 last[next[vi]] = last[vi];
1436 }
1437
1438 /* ------remove inactive items from element list of vi */
1439 L3:
1440 ls = vi;
1441 L4:
1442 s = ls;
1443 ls = l[s];
1444 if (ls == 0) {
1445 goto L6;
1446 }
1447 es = v[ls];
1448 if (mark[es] < tag) {
1449 goto L5;
1450 }
1451 free = ls;
1452 l[s] = l[ls];
1453 ls = s;
1454 L5:
1455 goto L4;
1456
1457 /* ------if vi is interior vertex, then remove from list and eliminate
1458 */
1459 L6:
1460 lvi = l[vi];
1461 if (lvi != 0) {
1462 goto L7;
1463 }
1464 l[i] = l[li];
1465 li = i;
1466
1467 ++(*k);
1468 next[vi] = -(*k);
1469 --last[*ek];
1470 goto L11;
1471
1472 /* ------else ... */
1473 /* --------classify vertex vi */
1474 L7:
1475 if (l[lvi] != 0) {
1476 goto L9;
1477 }
1478 evi = v[lvi];
1479 if (next[evi] >= 0) {
1480 goto L9;
1481 }
1482 if (mark[evi] < 0) {
1483 goto L8;
1484 }
1485
1486 /* ----------if vi is prototype vertex, then mark as such, initialize
1487 */
1488 /* ----------overlap count for corresponding element, and move vi to e
1489 nd */
1490 /* ----------of boundary list */
1491 last[vi] = evi;
1492 mark[evi] = -1;
1493 l[*tail] = li;
1494 *tail = li;
1495 l[i] = l[li];
1496 li = i;
1497 goto L10;
1498
1499 /* ----------else if vi is duplicate vertex, then mark as such and adj
1500 ust */
1501 /* ----------overlap count for corresponding element */
1502 L8:
1503 last[vi] = 0;
1504 --mark[evi];
1505 goto L10;
1506
1507 /* ----------else mark vi to compute degree */
1508 L9:
1509 last[vi] = -(*ek);
1510
1511 /* --------insert ek in element list of vi */
1512 L10:
1513 v[free] = *ek;
1514 l[free] = l[vi];
1515 l[vi] = free;
1516 L11:
1517 ;
1518 }
1519
1520 /* ----terminate boundary list */
1521 L12:
1522 l[*tail] = 0;
1523
1524 return 0;
1525 } /* mdp_ */
1526
mdm_(long * vk,long * tail,long * v,long * l,long * last,long * next,long * mark)1527 int mdm_(long *vk, long *tail, long *v, long *l,
1528 long *last, long *next, long *mark)
1529 {
1530 /* System generated locals */
1531 long i__1;
1532 static long equiv_0[1];
1533
1534 /* Local variables */
1535 long b, s, lb;
1536 #define es (equiv_0)
1537 long vb, ls;
1538 #define vs (equiv_0)
1539 long blpmax, tag, blp;
1540
1541 /* lll. optimize */
1542 /*
1543 */
1544 /* mdm -- form element from uneliminated neighbors of vk */
1545 /*
1546 */
1547
1548 /* ----initialize tag and list of uneliminated neighbors */
1549 /* Parameter adjustments */
1550 --mark;
1551 --next;
1552 --last;
1553 --l;
1554 --v;
1555
1556 /* Function Body */
1557 tag = mark[*vk];
1558 *tail = *vk;
1559
1560 /* ----for each vertex/element vs/es in element list of vk */
1561 ls = l[*vk];
1562 L1:
1563 s = ls;
1564 if (s == 0) {
1565 goto L5;
1566 }
1567 ls = l[s];
1568 *vs = v[s];
1569 if (next[*vs] < 0) {
1570 goto L2;
1571 }
1572
1573 /* ------if vs is uneliminated vertex, then tag and append to list of */
1574 /* ------uneliminated neighbors */
1575 mark[*vs] = tag;
1576 l[*tail] = s;
1577 *tail = s;
1578 goto L4;
1579
1580 /* ------if es is active element, then ... */
1581 /* --------for each vertex vb in boundary list of element es */
1582 L2:
1583 lb = l[*es];
1584 blpmax = last[*es];
1585 i__1 = blpmax;
1586 for (blp = 1; blp <= i__1; ++blp) {
1587 b = lb;
1588 lb = l[b];
1589 vb = v[b];
1590
1591 /* ----------if vb is untagged vertex, then tag and append to list of
1592 */
1593 /* ----------uneliminated neighbors */
1594 if (mark[vb] >= tag) {
1595 goto L3;
1596 }
1597 mark[vb] = tag;
1598 l[*tail] = b;
1599 *tail = b;
1600 L3:
1601 ;
1602 }
1603
1604 /* --------mark es inactive */
1605 mark[*es] = tag;
1606
1607 L4:
1608 goto L1;
1609
1610 /* ----terminate list of uneliminated neighbors */
1611 L5:
1612 l[*tail] = 0;
1613
1614 return 0;
1615 } /* mdm_ */
1616
1617 #undef vs
1618 #undef es
1619
1620
mdi_(long * n,long * ia,long * ja,long * max_,long * v,long * l,long * head,long * last,long * next,long * mark,long * tag,long * flag_)1621 int mdi_(long *n, long *ia, long *ja, long *max_,
1622 long *v, long *l, long *head, long *last,
1623 long *next, long *mark, long *tag, long *flag_)
1624 {
1625 /* System generated locals */
1626 long i__1, i__2, i__3;
1627
1628 /* Local variables */
1629 long jmin, jmax, kmax, j, k, vi, vj, nextvi, dvi, sfs, lvk;
1630
1631 /* lll. optimize */
1632 /*
1633 */
1634 /* mdi -- initialization */
1635 /*
1636 */
1637
1638 /* initialize degrees, element lists, and degree lists */
1639 /* Parameter adjustments */
1640 --mark;
1641 --next;
1642 --last;
1643 --head;
1644 --l;
1645 --v;
1646 --ja;
1647 --ia;
1648
1649 /* Function Body */
1650 i__1 = *n;
1651 for (vi = 1; vi <= i__1; ++vi) {
1652 mark[vi] = 1;
1653 l[vi] = 0;
1654 /* L1: */
1655 head[vi] = 0;
1656 }
1657 sfs = *n + 1;
1658
1659 /* create nonzero structure */
1660 /* for each nonzero entry a(vi,vj) */
1661 i__1 = *n;
1662 for (vi = 1; vi <= i__1; ++vi) {
1663 jmin = ia[vi];
1664 jmax = ia[vi + 1] - 1;
1665 if (jmin > jmax) {
1666 goto L6;
1667 }
1668 i__2 = jmax;
1669 for (j = jmin; j <= i__2; ++j) {
1670 vj = ja[j];
1671 if ((i__3 = vj - vi) < 0) {
1672 goto L2;
1673 } else if (i__3 == 0) {
1674 goto L5;
1675 } else {
1676 goto L4;
1677 }
1678
1679 /* ------if a(vi,vj) is in strict lower triangle */
1680 /* ------check for previous occurrence of a(vj,vi) */
1681 L2:
1682 lvk = vi;
1683 kmax = mark[vi] - 1;
1684 if (kmax == 0) {
1685 goto L4;
1686 }
1687 i__3 = kmax;
1688 for (k = 1; k <= i__3; ++k) {
1689 lvk = l[lvk];
1690 if (v[lvk] == vj) {
1691 goto L5;
1692 }
1693
1694 }
1695 /* ----for unentered entries a(vi,vj) */
1696 L4:
1697 if (sfs >= *max_) {
1698 goto L101;
1699 }
1700
1701 /* ------enter vj in element list for vi */
1702 ++mark[vi];
1703 v[sfs] = vj;
1704 l[sfs] = l[vi];
1705 l[vi] = sfs;
1706 ++sfs;
1707
1708 /* ------enter vi in element list for vj */
1709 ++mark[vj];
1710 v[sfs] = vi;
1711 l[sfs] = l[vj];
1712 l[vj] = sfs;
1713 ++sfs;
1714 L5:
1715 ;
1716 }
1717 L6:
1718 ;
1719 }
1720
1721 /* create degree lists and initialize mark vector */
1722 i__1 = *n;
1723 for (vi = 1; vi <= i__1; ++vi) {
1724 dvi = mark[vi];
1725 next[vi] = head[dvi];
1726 head[dvi] = vi;
1727 last[vi] = -dvi;
1728 nextvi = next[vi];
1729 if (nextvi > 0) {
1730 last[nextvi] = vi;
1731 }
1732 /* L7: */
1733 mark[vi] = *tag;
1734 }
1735
1736 return 0;
1737
1738 /* error- insufficient storage */
1739 L101:
1740 *flag_ = *n * 9 + vi;
1741 return 0;
1742 } /* mdi_ */
1743
1744
1745 /* -----------------------------------------------------------------------------
1746 jgroup_
1747
1748 this subroutine constructs groupings of the column indices of
1749 the jacobian matrix, used in the numerical evaluation of the
1750 jacobian by finite differences.
1751
1752 input..
1753 n = the order of the matrix.
1754 ia,ja = sparse structure descriptors of the matrix by rows.
1755 maxg = length of available storate in the igp array.
1756
1757 output..
1758 ngrp = number of groups.
1759 jgp = array of length n containing the column indices by groups.
1760 igp = pointer array of length ngrp + 1 to the locations in jgp
1761 of the beginning of each group.
1762 ier = error indicator. ier = 0 if no error occurred, or 1 if
1763 maxg was insufficient.
1764
1765 incl and jdone are working arrays of length n.
1766 */
1767
jgroup_(long * n,long * ia,long * ja,long * maxg,long * ngrp,long * igp,long * jgp,long * incl,long * jdone,long * ier)1768 int jgroup_(long *n, long *ia, long *ja, long *
1769 maxg, long *ngrp, long *igp, long *jgp,
1770 long *incl, long *jdone, long *ier)
1771 {
1772 /* Local variables */
1773 long ncol, kmin, kmax, i, j, k, ng;
1774
1775 /* Parameter adjustments */
1776 --jdone;
1777 --incl;
1778 --jgp;
1779 --igp;
1780 --ja;
1781 --ia;
1782
1783 /* Function Body */
1784 *ier = 0;
1785 for (j = 1; j <= *n; ++j) jdone[j] = 0;
1786
1787 ncol = 1;
1788 for (ng = 1; ng <= *maxg; ++ng) {
1789 igp[ng] = ncol;
1790 for (i = 1; i <= *n; ++i) incl[i] = 0;
1791
1792 for (j = 1; j <= *n; ++j) {
1793 /* reject column j if it is already in a group */
1794 if (jdone[j] == 1) goto L50;
1795
1796 kmin = ia[j];
1797 kmax = ia[j + 1] - 1;
1798 for (k = kmin; k <= kmax; ++k) {
1799 /* reject column j if it overlaps any column already in this group */
1800 i = ja[k];
1801 if (incl[i] == 1) goto L50;
1802 }
1803
1804 /* accept column j into group ng */
1805 jgp[ncol] = j;
1806 ++ncol;
1807 jdone[j] = 1;
1808 for (k = kmin; k <= kmax; ++k) {
1809 i = ja[k];
1810 incl[i] = 1;
1811 }
1812
1813 L50: ;
1814 }
1815 /* stop if this group is empty (grouping is complete) */
1816 if (ncol == igp[ng]) goto L70;
1817 }
1818
1819 /* error return if not all columns were chosen (maxg too small) */
1820 if (ncol <= *n) goto L80;
1821
1822 ng = *maxg;
1823
1824 L70:
1825 *ngrp = ng - 1;
1826 return 0;
1827
1828 L80:
1829 *ier = 1;
1830 return 0;
1831
1832 } /* jgroup_ */
1833
1834
1835 /* -----------------------------------------------------------------------------
1836 ewset_
1837
1838 this subroutine sets the error weight vector ewt according to
1839 ewt(i) = rtol(i)*abs(ycur(i)) + atol(i), i = 1,...,n,
1840 with the subscript on rtol and/or atol possibly replaced by 1 above,
1841 depending on the value of itol.
1842 */
1843
ewset_(long * n,long * itol,double * rtol,double * atol,double * ycur,double * ewt)1844 int ewset_(long *n, long *itol, double *rtol,
1845 double *atol, double *ycur, double *ewt)
1846 {
1847 long i;
1848
1849 /* Function Body */
1850 switch (*itol) {
1851 case 1:
1852 for (i = 0; i < *n; ++i)
1853 ewt[i] = rtol[0] * fabs(ycur[i]) + atol[0];
1854 break;
1855
1856 case 2:
1857 for (i = 0; i < *n; ++i)
1858 ewt[i] = rtol[0] * fabs(ycur[i]) + atol[i];
1859 break;
1860
1861 case 3:
1862 for (i = 0; i < *n; ++i)
1863 ewt[i] = rtol[i] * fabs(ycur[i]) + atol[0];
1864 break;
1865
1866 case 4:
1867 for (i = 0; i < *n; ++i)
1868 ewt[i] = rtol[i] * fabs(ycur[i]) + atol[i];
1869 break;
1870 }
1871
1872 return 0;
1873
1874 } /* ewset_ */
1875
1876 /* end */
1877