subroutine nsfc * (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju, * q, ira,jra, irac, irl,jrl, iru,jru, flag) clll. optimize c*** subroutine nsfc c*** symbolic ldu-factorization of nonsymmetric sparse matrix c (compressed pointer storage) c c c input variables.. n, r, ic, ia, ja, jlmax, jumax. c output variables.. il, jl, ijl, iu, ju, iju, flag. c c parameters used internally.. c nia - q - suppose m* is the result of reordering m. if c - processing of the ith row of m* (hence the ith c - row of u) is being done, q(j) is initially c - nonzero if m*(i,j) is nonzero (j.ge.i). since c - values need not be stored, each entry points to the c - next nonzero and q(n+1) points to the first. n+1 c - indicates the end of the list. for example, if n=9 c - and the 5th row of m* is c - 0 x x 0 x 0 0 x 0 c - then q will initially be c - a a a a 8 a a 10 5 (a - arbitrary). c - as the algorithm proceeds, other elements of q c - are inserted in the list because of fillin. c - q is used in an analogous manner to compute the c - ith column of l. c - size = n+1. c nia - ira, - vectors used to find the columns of m. at the kth c nia - jra, step of the factorization, irac(k) points to the c nia - irac head of a linked list in jra of row indices i c - such that i .ge. k and m(i,k) is nonzero. zero c - indicates the end of the list. ira(i) (i.ge.k) c - points to the smallest j such that j .ge. k and c - m(i,j) is nonzero. c - size of each = n. c nia - irl, - vectors used to find the rows of l. at the kth step c nia - jrl of the factorization, jrl(k) points to the head c - of a linked list in jrl of column indices j c - such j .lt. k and l(k,j) is nonzero. zero c - indicates the end of the list. irl(j) (j.lt.k) c - points to the smallest i such that i .ge. k and c - l(i,j) is nonzero. c - size of each = n. c nia - iru, - vectors used in a manner analogous to irl and jrl c nia - jru to find the columns of u. c - size of each = n. c c internal variables.. c jlptr - points to the last position used in jl. c juptr - points to the last position used in ju. c jmin,jmax - are the indices in a or u of the first and last c elements to be examined in a given row. c for example, jmin=ia(k), jmax=ia(k+1)-1. c integer cend, qm, rend, rk, vj integer ia(1), ja(1), ira(1), jra(1), il(1), jl(1), ijl(1) integer iu(1), ju(1), iju(1), irl(1), jrl(1), iru(1), jru(1) integer r(1), ic(1), q(1), irac(1), flag c c ****** initialize pointers **************************************** np1 = n + 1 jlmin = 1 jlptr = 0 il(1) = 1 jumin = 1 juptr = 0 iu(1) = 1 do 1 k=1,n irac(k) = 0 jra(k) = 0 jrl(k) = 0 1 jru(k) = 0 c ****** initialize column pointers for a *************************** do 2 k=1,n rk = r(k) iak = ia(rk) if (iak .ge. ia(rk+1)) go to 101 jaiak = ic(ja(iak)) if (jaiak .gt. k) go to 105 jra(k) = irac(jaiak) irac(jaiak) = k 2 ira(k) = iak c c ****** for each column of l and row of u ************************** do 41 k=1,n c c ****** initialize q for computing kth column of l ***************** q(np1) = np1 luk = -1 c ****** by filling in kth column of a ****************************** vj = irac(k) if (vj .eq. 0) go to 5 3 qm = np1 4 m = qm qm = q(m) if (qm .lt. vj) go to 4 if (qm .eq. vj) go to 102 luk = luk + 1 q(m) = vj q(vj) = qm vj = jra(vj) if (vj .ne. 0) go to 3 c ****** link through jru ******************************************* 5 lastid = 0 lasti = 0 ijl(k) = jlptr i = k 6 i = jru(i) if (i .eq. 0) go to 10 qm = np1 jmin = irl(i) jmax = ijl(i) + il(i+1) - il(i) - 1 long = jmax - jmin if (long .lt. 0) go to 6 jtmp = jl(jmin) if (jtmp .ne. k) long = long + 1 if (jtmp .eq. k) r(i) = -r(i) if (lastid .ge. long) go to 7 lasti = i lastid = long c ****** and merge the corresponding columns into the kth column **** 7 do 9 j=jmin,jmax vj = jl(j) 8 m = qm qm = q(m) if (qm .lt. vj) go to 8 if (qm .eq. vj) go to 9 luk = luk + 1 q(m) = vj q(vj) = qm qm = vj 9 continue go to 6 c ****** lasti is the longest column merged into the kth ************ c ****** see if it equals the entire kth column ********************* 10 qm = q(np1) if (qm .ne. k) go to 105 if (luk .eq. 0) go to 17 if (lastid .ne. luk) go to 11 c ****** if so, jl can be compressed ******************************** irll = irl(lasti) ijl(k) = irll + 1 if (jl(irll) .ne. k) ijl(k) = ijl(k) - 1 go to 17 c ****** if not, see if kth column can overlap the previous one ***** 11 if (jlmin .gt. jlptr) go to 15 qm = q(qm) do 12 j=jlmin,jlptr if (jl(j) - qm) 12, 13, 15 12 continue go to 15 13 ijl(k) = j do 14 i=j,jlptr if (jl(i) .ne. qm) go to 15 qm = q(qm) if (qm .gt. n) go to 17 14 continue jlptr = j - 1 c ****** move column indices from q to jl, update vectors *********** 15 jlmin = jlptr + 1 ijl(k) = jlmin if (luk .eq. 0) go to 17 jlptr = jlptr + luk if (jlptr .gt. jlmax) go to 103 qm = q(np1) do 16 j=jlmin,jlptr qm = q(qm) 16 jl(j) = qm 17 irl(k) = ijl(k) il(k+1) = il(k) + luk c c ****** initialize q for computing kth row of u ******************** q(np1) = np1 luk = -1 c ****** by filling in kth row of reordered a *********************** rk = r(k) jmin = ira(k) jmax = ia(rk+1) - 1 if (jmin .gt. jmax) go to 20 do 19 j=jmin,jmax vj = ic(ja(j)) qm = np1 18 m = qm qm = q(m) if (qm .lt. vj) go to 18 if (qm .eq. vj) go to 102 luk = luk + 1 q(m) = vj q(vj) = qm 19 continue c ****** link through jrl, ****************************************** 20 lastid = 0 lasti = 0 iju(k) = juptr i = k i1 = jrl(k) 21 i = i1 if (i .eq. 0) go to 26 i1 = jrl(i) qm = np1 jmin = iru(i) jmax = iju(i) + iu(i+1) - iu(i) - 1 long = jmax - jmin if (long .lt. 0) go to 21 jtmp = ju(jmin) if (jtmp .eq. k) go to 22 c ****** update irl and jrl, ***************************************** long = long + 1 cend = ijl(i) + il(i+1) - il(i) irl(i) = irl(i) + 1 if (irl(i) .ge. cend) go to 22 j = jl(irl(i)) jrl(i) = jrl(j) jrl(j) = i 22 if (lastid .ge. long) go to 23 lasti = i lastid = long c ****** and merge the corresponding rows into the kth row ********** 23 do 25 j=jmin,jmax vj = ju(j) 24 m = qm qm = q(m) if (qm .lt. vj) go to 24 if (qm .eq. vj) go to 25 luk = luk + 1 q(m) = vj q(vj) = qm qm = vj 25 continue go to 21 c ****** update jrl(k) and irl(k) *********************************** 26 if (il(k+1) .le. il(k)) go to 27 j = jl(irl(k)) jrl(k) = jrl(j) jrl(j) = k c ****** lasti is the longest row merged into the kth *************** c ****** see if it equals the entire kth row ************************ 27 qm = q(np1) if (qm .ne. k) go to 105 if (luk .eq. 0) go to 34 if (lastid .ne. luk) go to 28 c ****** if so, ju can be compressed ******************************** irul = iru(lasti) iju(k) = irul + 1 if (ju(irul) .ne. k) iju(k) = iju(k) - 1 go to 34 c ****** if not, see if kth row can overlap the previous one ******** 28 if (jumin .gt. juptr) go to 32 qm = q(qm) do 29 j=jumin,juptr if (ju(j) - qm) 29, 30, 32 29 continue go to 32 30 iju(k) = j do 31 i=j,juptr if (ju(i) .ne. qm) go to 32 qm = q(qm) if (qm .gt. n) go to 34 31 continue juptr = j - 1 c ****** move row indices from q to ju, update vectors ************** 32 jumin = juptr + 1 iju(k) = jumin if (luk .eq. 0) go to 34 juptr = juptr + luk if (juptr .gt. jumax) go to 106 qm = q(np1) do 33 j=jumin,juptr qm = q(qm) 33 ju(j) = qm 34 iru(k) = iju(k) iu(k+1) = iu(k) + luk c c ****** update iru, jru ******************************************** i = k 35 i1 = jru(i) if (r(i) .lt. 0) go to 36 rend = iju(i) + iu(i+1) - iu(i) if (iru(i) .ge. rend) go to 37 j = ju(iru(i)) jru(i) = jru(j) jru(j) = i go to 37 36 r(i) = -r(i) 37 i = i1 if (i .eq. 0) go to 38 iru(i) = iru(i) + 1 go to 35 c c ****** update ira, jra, irac ************************************** 38 i = irac(k) if (i .eq. 0) go to 41 39 i1 = jra(i) ira(i) = ira(i) + 1 if (ira(i) .ge. ia(r(i)+1)) go to 40 irai = ira(i) jairai = ic(ja(irai)) if (jairai .gt. i) go to 40 jra(i) = irac(jairai) irac(jairai) = i 40 i = i1 if (i .ne. 0) go to 39 41 continue c ijl(n) = jlptr iju(n) = juptr flag = 0 return c c ** error.. null row in a 101 flag = n + rk return c ** error.. duplicate entry in a 102 flag = 2*n + rk return c ** error.. insufficient storage for jl 103 flag = 3*n + k return c ** error.. null pivot 105 flag = 5*n + k return c ** error.. insufficient storage for ju 106 flag = 6*n + k return end