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