subroutine rscf3 (lftblk, array, right, prodx1, prodx2, nrwblk, * nbloks, rgtblk, pivot, iflag, nparts, blaws) c double precision lftblk(nrwblk,1), array(nrwblk,2*nrwblk,1), * right(nrwblk,nrwblk,1), prodx1(nrwblk,nrwblk,1), * prodx2(nrwblk,1), rgtblk(nrwblk,1), blaws(nrwblk,1) integer nrwblk, nbloks, pivot(1), iflag, nparts c ***----------------------------------------------------------*** c * The following notation is used in the comments: * c * * c * B_a, B_b <=> lftblk, rgtblk * c * V_k, W_k <=> array(1,1,k), array(1,nrwblk+1,k) * c * R_k <=> right(1,1,k) * c * S_k <=> prodx1(1,1,k) * c * T <=> prodx2(1,1) * c * * c * In addition, the affix '/'' designates that the * c * matrix was (or at least could have been) transformed * c * at the first/second level of the factorization. * c * (In order to be consistent with other variants of * c * this algorithm, all transformations to B_a and B_b * c * occur at the third level of the factorization.) * c ***----------------------------------------------------------*** integer nsquar, info c nsquar = nrwblk**2 c ***----------------------------------------------------------*** c * R_nbloks <- LUfact(W_nbloks) * c * * c * [ Note: The LU factorization of W_nbloks is stored * c * in R_nbloks for use below and for processing * c * subsequent right hand sides. The pivot indices * c * are stored in pivot(1..nrwblk). ] * c ***----------------------------------------------------------*** call DCOPY(nsquar,array(1,nrwblk+1,nbloks),1, * right(1,1,nbloks),1) call DGETRF(nrwblk,nrwblk,right(1,1,nbloks),nrwblk, * pivot,iflag) if (iflag .ne. 0) return c ***----------------------------------------------------------*** c * The content of part of the third-level block-array * c * depends on whether or not paritioning was done. * c ***----------------------------------------------------------*** if (nparts .gt. 1) then c ***----------------------------------------------------------*** c * If there is more than one partition, the third-level * c * block-array is of the form [B_a'' B_b * c * T'' W_nbloks], * c * where B_a'' = B_a' + B_a' S_1'' * c * = B_a' (I + S_1'') * c * = (B_a + B_a V_1') (I + S_1'') * c * = B_a (I + V_1') (I + S_1'') * c * = B_a (I + V_1' + S_1'' + V_1' S_1''), * c * * c * and T'' = +/-S_nparts' S_nparts-1''... S_1''. * c * * c * Instead of directly factoring this 2x2 block array, * c * the reduced array [B_a'' - B_b W_nbloks^-1 T''] is * c * factored and its factorization is stored at the base * c * of the last partition in W_nbloks. * c * * c * [ Note: This cannot be done in the slf_lu or slf_qr * c * algorithms since W_nbloks is often numerically * c * singular at this level of the factorization. ] * c ***----------------------------------------------------------*** c c ***----------------------------------------------------------*** c * W_nbloks <- B_a (I + V_1' + S_1'' + V_1' S_1'') * c ***----------------------------------------------------------*** call DCOPY(nsquar,prodx1,1,blaws,1) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * 1.d0,array,nrwblk,prodx1,nrwblk, * 1.d0,blaws,nrwblk) call DAXPY(nsquar,1.d0,array,1,blaws,1) call maddi('+',nrwblk,blaws) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * 1.d0,lftblk,nrwblk,blaws,nrwblk, * 0.d0,array(1,nrwblk+1,nbloks),nrwblk) c ***----------------------------------------------------------*** c * W_nbloks <- (B_a'' - B_b W_nbloks^-1 T'') * c ***----------------------------------------------------------*** call DCOPY(nsquar,prodx2,1,blaws,1) call DGETRS('N',nrwblk,nrwblk,right(1,1,nbloks), * nrwblk,pivot,blaws,nrwblk,info) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * -1.d0,rgtblk,nrwblk,blaws,nrwblk, * 1.d0,array(1,nrwblk+1,nbloks),nrwblk) c ***----------------------------------------------------------*** c * W_nbloks <- LUfact(B_a'' - B_b W_nbloks^-1 T'') * c ***----------------------------------------------------------*** call DGETRF(nrwblk,nrwblk,array(1,nrwblk+1,nbloks),nrwblk, * pivot(nbloks*nrwblk+1),iflag) else c ***----------------------------------------------------------*** c * If there is only one partition, the third-level * c * block-array is of the form [B_a' B_b * c * S_1' W_nbloks], * c * where B_a' = B_a + B_a V_1', * c * * c * and S_1' = +/-V_nbloks V_nbloks-1'... V_1'. * c * * c * Instead of directly factoring this 2x2 block array, * c * the reduced array [B_a' - B_b W_nbloks^-1 S_1'] is * c * factored and its factorization is stored in W_nbloks. * c ***----------------------------------------------------------*** c c ***----------------------------------------------------------*** c * W_nbloks <- (B_a + B_a V_1') * c ***----------------------------------------------------------*** call DCOPY(nsquar,lftblk,1,array(1,nrwblk+1,nbloks),1) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * 1.d0,lftblk,nrwblk,array,nrwblk, * 1.d0,array(1,nrwblk+1,nbloks),nrwblk) c ***----------------------------------------------------------*** c * W_nbloks <- (B_a' - B_b W_nbloks^-1 S_1') * c ***----------------------------------------------------------*** call DCOPY(nsquar,prodx1,1,blaws,1) call DGETRS('N',nrwblk,nrwblk,right(1,1,nbloks), * nrwblk,pivot,blaws,nrwblk,info) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * -1.d0,rgtblk,nrwblk,blaws,nrwblk, * 1.d0,array(1,nrwblk+1,nbloks),nrwblk) c ***----------------------------------------------------------*** c * W_nbloks <- LUfact(B_a' - B_b W_nbloks^-1 S_1') * c ***----------------------------------------------------------*** call DGETRF(nrwblk,nrwblk,array(1,nrwblk+1,nbloks),nrwblk, * pivot(nbloks*nrwblk+1),iflag) endif return end