1      subroutine y12mcf(n,z,a,snr,nn,rnr,nn1,pivot,b,ha,iha, aflag,iflag
2     *,ifail)
3c
4c  systens of linear equations are solved by use of sparse matrix tech-
5c  nique and by gaussian elimination.
6c
7      implicit double precision(a-b,g,p,t-y),integer(c,f,h-n,r-s,z)
8      double precision a(nn),b(n),pivot(n),aflag(8)
9c
10c  information which is necessary to begin the elimination is stored.
11c
12      integer snr(nn),rnr(nn1),ha(iha,11), iflag(10)
13      ifail=0
14      if(iflag(1).ne.-1)ifail=2
15      if(aflag(1).lt.1.0d0)aflag(1)=1.0005 d0
16      if(aflag(3).lt.1.0d+5)aflag(3)=1.0d+5
17      if(aflag(4).lt.0.0d0)aflag(4)=-aflag(4)
18      if(iflag(2).lt.1)ifail=19
19      if(iflag(3).lt.0.or.iflag(3).gt.2)ifail=20
20      if(iflag(5).lt.1.or.iflag(5).gt.3)ifail=21
21      if(iflag(5).eq.3)ifail=22
22      if(ifail.gt.0)go to 1110
23      snr(z+1)=0
24      rnr(z+1)=0
25      n8=n+1
26      n7=n-1
27      u=aflag(1)
28      grmin=aflag(4)*aflag(6)
29c
30c  use the information about fill-ins if it is possible.
31c
32      zz=z
33      nr=n*n
34      if(iflag(4).ne.2)go to 100
35      if(iflag(10).gt.nn)go to 50
36      l1=iflag(10)
37      l5=l1+1
38      if(l5.le.nn)snr(l5)=0
39      do 40 i=1,n
40      l=n8-i
41      l2=ha(l,3)+1
42      l3=l2-ha(l,1)
43      do 10 j=1,l3
44      snr(l5-j)=snr(l2-j)
45   10 a(l5-j)=a(l2-j)
46      ha(l,3)=l1
47      ha(l,1)=l5-l3
48      l6=l1-l3
49      l5=l5-ha(l,9)
50      if(l5.gt.l6)go to 30
51      do 20 j=l5,l6
52   20 snr(j)=0
53   30 continue
54   40 l1=l5-1
55   50 if(iflag(9).gt.nn1)go to 100
56      l2=iflag(9)
57      l5=l2+1
58      if(l5.le.nn1)rnr(l5)=0
59      do 90 i=1,n
60      l=n8-i
61      l1=ha(l,6)+1
62      l4=l1-ha(l,4)
63      do 60 j=1,l4
64   60 rnr(l5-j)=rnr(l1-j)
65      ha(l,4)=l5-l4
66      ha(l,6)=l2
67      l6=l2-l4
68      l5=l5-ha(l,10)
69      if(l5.gt.l6)go to 80
70      do 70 j=l5,l6
71   70 rnr(j)=0
72   80 continue
73   90 l2=l5-1
74  100 r4=ha(n,3)
75      r5=ha(n,6)
76      aflag(7)=aflag(6)
77      aflag(8)=aflag(6)
78      do 110 i=1,n
79      pivot(i)=0.0 d0
80      ha(i,2)=ha(i,1)
81  110 ha(i,5)=ha(i,4)
82      index=ha(n,8)
83c
84c  start of gaussian elimination.
85c
86      slut=ha(index,3)-ha(index,2)+1
87      do 950 i=1,n7
88      rr3=ha(i,2)
89      rr4=ha(i,3)
90      c1=ha(i,4)
91      cr4=ha(i,6)
92      if(iflag(3).eq.0)go to 350
93      if(iflag(4).ne.2)go to 120
94      rrow=ha(i,7)
95      rcoll=ha(i,8)
96      go to 220
97  120 l4=ha(i,8)
98      if(iflag(3).eq.1)go to 130
99      rrow=l4
100      rcoll=rrow
101      rpivot=i
102      go to 170
103  130 r=nr
104      v=0.0 d0
105      index=iflag(2)
106      do 160 kk=1,index
107      l1=i-1+kk
108      if(l1.gt.n)go to 170
109      j=ha(l1,8)
110      r7=ha(j,2)
111      r8=ha(j,3)
112      r9=r8-r7
113      t=0.0 d0
114      do 140 k=r7,r8
115      td=dabs(a(k))
116  140 if(t.lt.td)t=td
117      t=t/u
118      do 160 k=r7,r8
119      td=dabs(a(k))
120      if(td.lt.t)go to 150
121      r6=snr(k)
122      r3=r9*(ha(r6,6)-ha(r6,5))
123      if(r3.gt.r)go to 150
124      if(r3.lt.r)go to 151
125      if(v.ge.td)go to 150
126  151 v=td
127      rrow=j
128      rcoll=r6
129      r=r3
130      rpivot=l1
131  150 continue
132  160 continue
133  170 r3=ha(rcoll,10)
134      ha(rcoll,10)=ha(i,10)
135      ha(i,10)=r3
136      r3=ha(rrow,9)
137      ha(rrow,9)=ha(i,9)
138c
139c  remove the pivot row of the list where the rows are ordered by
140c  increasing numbers of non-zero elements.
141c
142      ha(i,9)=r3
143      l1=0
144      l=i
145      l2=ha(l4,3)-ha(l4,2)+1
146  180 l=l+1
147      if(l2.gt.l1)ha(l2,11)=l
148      if(l.gt.n)go to 190
149      l5=ha(l,8)
150      l3=ha(l5,3)-ha(l5,2)+1
151      if(rpivot.lt.l)go to 190
152      ha(l4,7)=l
153      ha(l,8)=l4
154      l4=l5
155      l1=l2
156      l2=l3
157      l3=n8
158      go to 180
159  190 if(l2.eq.l1)go to 200
160      if(l3.eq.l2)go to 200
161      ha(l2,11)=0
162  200 l5=ha(i,7)
163      if(rrow.eq.i)go to 210
164      ha(l5,8)=rrow
165      ha(rrow,7)=l5
166  210 ha(i,7)=rrow
167c
168c  row interchanges.
169c
170      ha(i,8)=rcoll
171  220 if(rrow.eq.i)go to 290
172      t=b(rrow)
173      b(rrow)=b(i)
174      b(i)=t
175      do 250 j=rr3,rr4
176      l1=snr(j)
177      r=ha(l1,5)-1
178      r10=ha(l1,6)
179  240 r=r+1
180      if(rnr(r).ne.i)go to 240
181      rnr(r)=rnr(r10)
182  250 rnr(r10)=rrow
183      rr3=ha(rrow,2)
184      rr4=ha(rrow,3)
185      do 270 j=rr3,rr4
186      l1=snr(j)
187      r=ha(l1,5)-1
188  260 r=r+1
189      if(rnr(r).ne.rrow)go to 260
190  270 rnr(r)=i
191      do 280 j=1,3
192      r3=ha(rrow,j)
193      ha(rrow,j)=ha(i,j)
194c
195c  column interchanges.
196c
197  280 ha(i,j)=r3
198  290 if(rcoll.eq.i)go to 350
199      do 310 j=c1,cr4
200      l1=rnr(j)
201      r=ha(l1,2)-1
202      r10=ha(l1,3)
203  300 r=r+1
204      if(snr(r).ne.i)go to 300
205      t=a(r10)
206      a(r10)=a(r)
207      a(r)=t
208      snr(r)=snr(r10)
209  310 snr(r10)=rcoll
210      c1=ha(rcoll,4)
211      cr4=ha(rcoll,6)
212      do 330 j=c1,cr4
213      l1=rnr(j)
214      r=ha(l1,2)-1
215  320 r=r+1
216      if(snr(r).ne.rcoll)go to 320
217  330 snr(r)=i
218      do 340 j=4,6
219      r3=ha(rcoll,j)
220      ha(rcoll,j)=ha(i,j)
221c
222c end of the interchanges.
223c the row ordered list and the column ordered list are prepared to
224c begin step i of the elimination.
225c
226  340 ha(i,j)=r3
227  350 r9=rr4-rr3
228      do 360 rr=rr3,rr4
229      if(snr(rr).eq.i)go to 370
230  360 continue
231      ifail=9
232      go to 1110
233  370 v=a(rr)
234      pivot(i)=v
235      td=dabs(v)
236      if(td.lt.aflag(8))aflag(8)=td
237      if(td.ge.grmin)go to 380
238      ifail=3
239      go to 1110
240  380 r2=ha(i,1)
241      a(rr)=a(rr3)
242      snr(rr)=snr(rr3)
243      a(rr3)=a(r2)
244      snr(rr3)=snr(r2)
245      snr(r2)=0
246      z=z-1
247      rr3=rr3+1
248      ha(i,2)=rr3
249      ha(i,1)=r2+1
250      cr3=ha(i,5)
251      if(r9.le.0)go to 431
252      do 430 j=rr3,rr4
253      index=snr(j)
254  430 pivot(index)=a(j)
255  431 r7=cr4-cr3+1
256      do 880 k=1,r7
257      r1=rnr(cr3-1+k)
258      if(r1.eq.i)go to 870
259      i1=ha(r1,1)
260      rr1=ha(r1,2)
261      rr2=ha(r1,3)
262      l2=rr2-rr1+1
263      l=rr1-1
264  390 l=l+1
265      if(snr(l).ne.i)go to 390
266      t=a(l)/v
267      if(iflag(5).eq.2)go to 400
268      a(l)=a(i1)
269      snr(l)=snr(i1)
270      snr(i1)=0
271      i1=i1+1
272      ha(r1,1)=i1
273      z=z-1
274      go to 410
275  400 a(l)=a(rr1)
276      a(rr1)=t
277      r3=snr(rr1)
278      snr(rr1)=snr(l)
279      snr(l)=r3
280  410 rr1=rr1+1
281      ha(r1,2)=rr1
282      b(r1)=b(r1)-b(i)*t
283      if(r9.le.0)go to 669
284      r=rr1
285      if(r.gt.rr2)go to 470
286      do 460 l=r,rr2
287      l1=snr(l)
288      td=pivot(l1)
289      if(td.eq.0.0d0)go to 450
290      pivot(l1)=0.0 d0
291      td=a(l)-td*t
292      a(l)=td
293      td1=dabs(td)
294      if(td1.gt.aflag(7))aflag(7)=td1
295c
296c  too small element is created.remove it from the lists.
297c
298      if(td1.gt.aflag(2))go to 450
299      z=z-1
300      a(l)=a(rr1)
301      snr(l)=snr(rr1)
302      a(rr1)=a(i1)
303      snr(rr1)=snr(i1)
304      snr(i1)=0
305      rr1=rr1+1
306      i1=i1+1
307      ha(r1,2)=rr1
308      ha(r1,1)=i1
309      r3=ha(l1,5)
310      r2=r3-1
311      l4=ha(l1,4)
312      l5=rnr(l4)
313      l6=rnr(r3)
314  440 r2=r2+1
315      if(rnr(r2).ne.r1)go to 440
316      rnr(r2)=l6
317      rnr(r3)=l5
318      rnr(l4)=0
319      ha(l1,5)=r3+1
320      ha(l1,4)=l4+1
321  450 continue
322  460 continue
323  470 continue
324      do 750 j=1,r9
325      r=rr3-1+j
326      r2=snr(r)
327      tol2=pivot(r2)
328      pivot(r2)=a(r)
329      if(tol2.eq.0.0d0)go to 740
330      tol3=-tol2*t
331      tol1=dabs(tol3)
332      if(tol1.lt.aflag(2))go to 740
333      c2=ha(r2,4)
334      cr2=ha(r2,6)
335      cr1=ha(r2,5)
336      lfr=rr2-i1+2
337      lfc=cr2-c2+2
338      if(iflag(4).ne.1)go to 480
339      if(lfr.gt.ha(r1,9))ha(r1,9)=lfr
340      if(lfc.gt.ha(r2,10))ha(r2,10)=lfc
341  480 if(i1.eq.1)go to 490
342      if(snr(i1-1).eq.0)go to 600
343  490 if(rr2.eq.nn)go to 500
344      if(snr(rr2+1).eq.0)go to 580
345c
346c  collection in row ordered list.
347c
348  500 r10=nn-lfr
349      if(r10.ge.r4)go to 560
350      iflag(6)=iflag(6)+1
351      do 520 jj=1,n
352      l1=ha(jj,3)
353      if(l1.lt.ha(jj,1))go to 510
354      ha(jj,3)=snr(l1)
355      snr(l1)=-jj
356  510 continue
357  520 continue
358      l3=0
359      l4=1
360      do 550 jj=1,r4
361      if(snr(jj).eq.0)go to 540
362      l3=l3+1
363      if(snr(jj).gt.0)go to 530
364      l5=-snr(jj)
365      snr(jj)=ha(l5,3)
366      ha(l5,3)=l3
367      l6=l4+ha(l5,2)-ha(l5,1)
368      ha(l5,2)=l6
369      ha(l5,1)=l4
370      l4=l3+1
371  530 a(l3)=a(jj)
372      snr(l3)=snr(jj)
373  540 continue
374  550 continue
375      r4=l3
376      snr(l3+1)=0
377      rr3=ha(i,2)
378      rr4=ha(i,3)
379      i1=ha(r1,1)
380      rr1=ha(r1,2)
381      r=rr3-1+j
382      if(r10.ge.r4)go to 560
383      ifail=5
384c
385c fill-in takes place in the row ordered list.
386c
387      go to 1110
388  560 r8=lfr-1
389      rr2=r4+lfr
390      if(r8.le.0)go to 579
391      l3=i1-1
392      do 570 ll=1,r8
393      l4=r4+ll
394      l5=l3+ll
395      a(l4)=a(l5)
396      snr(l4)=snr(l5)
397  570 snr(l5)=0
398  579 rr1=r4+rr1-i1+1
399      ha(r1,3)=rr2
400      ha(r1,2)=rr1
401      i1=r4+1
402      ha(r1,1)=i1
403      l1=rr2
404      go to 590
405  580 rr2=rr2+1
406      ha(r1,3)=rr2
407      l1=rr2
408      if(rr2.le.r4)go to 610
409  590 r4=rr2
410      if(r4.lt.nn)snr(r4+1)=0
411      go to 610
412  600 rr1=rr1-1
413      i1=i1-1
414      ha(r1,1)=i1
415      ha(r1,2)=rr1
416      l1=rr1
417      snr(i1)=snr(l1)
418      a(i1)=a(l1)
419  610 a(l1)=tol3
420      snr(l1)=snr(r)
421      td=dabs(a(l1))
422      if(td.gt.aflag(7))aflag(7)=td
423      z=z+1
424      if(iflag(8).lt.z) iflag(8)=z
425      if(c2.eq.1)go to 620
426      if(rnr(c2-1).eq.0)go to 720
427  620 if(cr2.eq.nn1)go to 630
428      if(rnr(cr2+1).eq.0)go to 700
429c
430c  collection in column ordered list.
431c
432  630 r10=nn1-lfc
433      if(r10.ge.r5)go to 680
434      iflag(7)=iflag(7)+1
435      do 640 jj=i,n
436      l1=ha(jj,6)
437      ha(jj,6)=rnr(l1)
438  640 rnr(l1)=-jj
439      l3=0
440      l4=1
441      do 670 jj=1,r5
442      if(rnr(jj).eq.0)go to 660
443      l3=l3+1
444      if(rnr(jj).gt.0)go to 650
445      l5=-rnr(jj)
446      rnr(jj)=ha(l5,6)
447      ha(l5,6)=l3
448      l6=l4+ha(l5,5)-ha(l5,4)
449      ha(l5,5)=l6
450      ha(l5,4)=l4
451      l4=l3+1
452  650 rnr(l3)=rnr(jj)
453  660 continue
454  670 continue
455      r5=l3
456      rnr(r5+1)=0
457      c2=ha(r2,4)
458      cr3=ha(i,5)
459      cr4=ha(i,6)
460      cr1=ha(r2,5)
461      if(r10.ge.r5)go to 680
462      ifail=6
463c
464c fill-in takes place in the column ordered list.
465c
466      go to 1110
467  680 r8=lfc-1
468      cr2=r5+lfc
469      if(r8.le.0)go to 699
470      l3=c2-1
471      do 690 l=1,r8
472      l4=r5+l
473      l5=l3+l
474      rnr(l4)=rnr(l5)
475  690 rnr(l5)=0
476  699 cr1=r5+cr1-c2+1
477      c2=r5+1
478      ha(r2,6)=cr2
479      ha(r2,4)=c2
480      ha(r2,5)=cr1
481      r=cr2
482      go to 710
483  700 cr2=cr2+1
484      ha(r2,6)=cr2
485      r=cr2
486      if(cr2.le.r5)go to 730
487  710 r5=cr2
488      if(r5.lt.nn1)rnr(r5+1)=0
489      go to 730
490  720 cr1=cr1-1
491      c2=c2-1
492      ha(r2,4)=c2
493      ha(r2,5)=cr1
494      r=cr1
495      rnr(c2)=rnr(r)
496  730 rnr(r)=r1
497  740 continue
498  750 continue
499  669 if(rr1.le.rr2)go to 760
500      ifail=7
501c
502c  update the information in the list where the rows are ordered by
503c  increasing numbers of the non-zero elements.
504c
505      go to 1110
506  760 if(iflag(4).eq.2)go to 870
507      if(iflag(3).eq.0)go to 870
508      l1=rr2-rr1+1
509      if(l1.eq.l2)go to 870
510      l6=ha(r1,7)
511      l4=ha(l2,11)
512      if(l1.gt.l2)go to 820
513      if(l6.gt.l4)go to 780
514      if(l4.eq.n)go to 770
515      l=ha(l4+1,8)
516      l5=ha(l,3)-ha(l,2)+1
517      if(l5.eq.l2)go to 790
518  770 ha(l2,11)=0
519      go to 800
520  780 l5=ha(l4,8)
521      l3=ha(l6,8)
522      ha(l4,8)=l3
523      ha(l6,8)=l5
524      ha(l5,7)=l6
525      ha(l3,7)=l4
526      l6=l4
527  790 ha(l2,11)=l4+1
528  800 if(l4.eq.i+1)go to 810
529      l=ha(l6-1,8)
530      l2=ha(l,3)-ha(l,2)+1
531      l4=ha(l2,11)
532      if(l1.lt.l2)go to 780
533  810 if(l1.ne.l2)ha(l1,11)=l6
534      go to 870
535  820 if(l6.gt.l4)go to 840
536      if(l4.eq.n)go to 830
537      l=ha(l4+1,8)
538      l5=ha(l,3)-ha(l,2)+1
539      if(l5.eq.l2)go to 840
540  830 ha(l2,11)=0
541  840 l2=l2+1
542      if(l2.le.slut)go to 850
543      l3=n
544      slut=l1
545      l2=l1
546      go to 860
547  850 l3=ha(l2,11)-1
548      if(l3.eq.-1)go to 840
549      if(l2.gt.l1)l2=l1
550  860 ha(l2,11)=l3
551      l4=ha(l3,8)
552      l7=ha(l6,8)
553      ha(l3,8)=l7
554      ha(l6,8)=l4
555      ha(l7,7)=l3
556      ha(l4,7)=l6
557      l6=l3
558      if(l2.lt.l1)go to 840
559  870 continue
560  880 continue
561      if(r9.le.0)go to 882
562      do 881 j=rr3,rr4
563      index=snr(j)
564  881 pivot(index)=0.0 d0
565  882 continue
566      cr3=ha(i,4)
567      do 890 j=cr3,cr4
568  890 rnr(j)=0
569      if(r9.le.0)go to 930
570      l2=ha(i,2)-1
571      do 920 ll=1,r9
572      r=snr(l2+ll)
573      r1=ha(r,5)
574      r2=ha(r,6)
575      if(r2.gt.r1)go to 900
576      ifail=8
577      go to 1110
578  900 ha(r,5)=r1+1
579      r3=r1-1
580  910 r3=r3+1
581      if(rnr(r3).ne.i)go to 910
582      rnr(r3)=rnr(r1)
583  920 rnr(r1)=i
584  930 aflag(5)=aflag(7)/aflag(6)
585      if(aflag(5).lt.aflag(3))go to 940
586      ifail=4
587      go to 1110
588  940 continue
589c
590c  preparation to begin the back substitution.
591c
592  950 continue
593      index=ha(n,2)
594      pivot(n)=a(index)
595      a(index)=0.0 d0
596      td=dabs(pivot(n))
597      if(td.gt.aflag(7))aflag(7)=td
598      if(td.lt.aflag(8))aflag(8)=td
599      if(td.gt.grmin)go to 960
600      ifail=3
601      go to 1110
602  960 if(iflag(4).ne.1)go to 1060
603      iflag(10)=ha(n,9)
604      iflag(9)=ha(n,10)
605      do 990 i=1,n7
606      r1=n-i
607      iflag(10)=iflag(10)+ha(r1,9)
608      iflag(9)=iflag(9)+ha(r1,10)
609      if(iflag(3).eq.0)go to 980
610      do 970 j=9,10
611      r2=ha(r1,j-2)
612      r6=ha(r2,j)
613      ha(r2,j)=ha(r1,j)
614  970 ha(r1,j)=r6
615  980 continue
616  990 continue
6171060  continue
618      aflag(5)=aflag(7)/aflag(6)
619      iflag(1)=-2
620 1110 z=zz
621      return
622      end
623