subroutine rscf2 (array, right, prodx1, prodx2, nrwblk, * pivot, iflag, minblk, remblk, nparts, blaws) c double precision array(nrwblk,2*nrwblk,1), * right(nrwblk,nrwblk,1), prodx1(nrwblk,nrwblk,1), * prodx2(nrwblk,1), blaws(nrwblk,1) integer nrwblk, pivot(1), iflag, minblk, remblk, nparts c ***----------------------------------------------------------*** c * The following notation is used in the comments: * c * * 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/is transformed at the first/second level * c * of the factorization. * c ***----------------------------------------------------------*** integer nsquar, kpart, nparts1, base, top, info c ***----------------------------------------------------------*** c * If there is only one partition, nothing needs to be * c * done at the second level of the factorization. * c ***----------------------------------------------------------*** if (nparts .eq. 1) return c ***----------------------------------------------------------*** c * Rescaling starts at the second-last block-row of * c * the second-level array. * c ***----------------------------------------------------------*** nsquar = nrwblk**2 nparts1 = nparts - 1 call partx(minblk,remblk,nparts1,base,top) c ***----------------------------------------------------------*** c * W_base'' <- (W_base' - S_nparts-1') * c ***----------------------------------------------------------*** call DAXPY(nsquar,-1.d0,prodx1(1,1,nparts1),1, * array(1,nrwblk+1,base),1) c ***----------------------------------------------------------*** c * W_base'' <- LUfact(W_base' - S_nparts-1') * c ***----------------------------------------------------------*** call DGETRF(nrwblk,nrwblk,array(1,nrwblk+1,base),nrwblk, * pivot(base*nrwblk+1),iflag) if (iflag .ne. 0) return c ***----------------------------------------------------------*** c * S_nparts-1'' <- (W_base' - S_nparts-1')^-1 S_nparts-1' * c ***----------------------------------------------------------*** call DGETRS('N',nrwblk,nrwblk,array(1,nrwblk+1,base), * nrwblk,pivot(base*nrwblk+1), * prodx1(1,1,nparts1),nrwblk,info) c ***----------------------------------------------------------*** c * T'' <- (S_nparts' S_nparts-1'')^T * c * * c * [ Notes: 1. The transpose of the product is accumulated * c * since Lapack's DGEMM is faster multiplying in * c * this mode when both matrices are dense. * c * * c * 2. After the last call to DGEMM, the resulting * c * matrix must be transposed once. ] * c ***----------------------------------------------------------*** call DGEMM('T','T',nrwblk,nrwblk,nrwblk, * 1.d0,prodx1(1,1,nparts1),nrwblk, * prodx1(1,1,nparts),nrwblk, * 0.d0,prodx2,nrwblk) c ***----------------------------------------------------------*** c * The second-level array is now processed sequentially * c * from the third-last block row to the top. * c ***----------------------------------------------------------*** do 10 kpart = nparts-2, 1, -1 call partx(minblk,remblk,kpart,base,top) c ***----------------------------------------------------------*** c * W_base'' <- W_base'(I + S_kpart+1'') - S_kpart' * c * * c * [ Note: A copy of W_base' was stored in R_base-1 during * c * the first-level factorization. ] * c ***----------------------------------------------------------*** call DCOPY(nsquar,prodx1(1,1,kpart+1),1,blaws,1) call maddi('+',nrwblk,blaws) call DCOPY(nsquar,prodx1(1,1,kpart),1, * array(1,nrwblk+1,base),1) call DGEMM('N','N',nrwblk,nrwblk,nrwblk, * 1.d0,right(1,1,base-1),nrwblk,blaws,nrwblk, * -1.d0,array(1,nrwblk+1,base),nrwblk) c ***----------------------------------------------------------*** c * W_base'' <- LUfact(W_base'(I + S_kpart+1'') - S_kpart') * c ***----------------------------------------------------------*** call DGETRF(nrwblk,nrwblk,array(1,nrwblk+1,base),nrwblk, * pivot(base*nrwblk+1),iflag) if (iflag .ne. 0) return c ***----------------------------------------------------------*** c * S_kpart'' <- (W_base'(I + S_kpart+1'') * c * - S_kpart')^-1 S_kpart' * c ***----------------------------------------------------------*** call DGETRS('N',nrwblk,nrwblk,array(1,nrwblk+1,base), * nrwblk,pivot(base*nrwblk+1), * prodx1(1,1,kpart),nrwblk,info) c ***----------------------------------------------------------*** c * T'' <- S_kpart''^T T'' * c ***----------------------------------------------------------*** call DGEMM('T','N',nrwblk,nrwblk,nrwblk, * 1.d0,prodx1(1,1,kpart),nrwblk, * prodx2,nrwblk,0.d0,blaws,nrwblk) call DCOPY(nsquar,blaws,1,prodx2,1) 10 continue c ***----------------------------------------------------------*** c * The final product must be transposed. * c ***----------------------------------------------------------*** call mtran(nrwblk,prodx2) c ***----------------------------------------------------------*** c * If there was an odd number of multiplications, * c * the final product also must be negated. * c ***----------------------------------------------------------*** if (mod(nparts1,2) .ne. 0) then call mnegv(nrwblk,prodx2) end if return end