! File: gsfct.f ! Date: 12/31/01 (downloaded fall 2001) ! Author: C. Paciorek ! Purpose: Cholesky factorization of positive definite matrices stored in sparse format; double precision version ! Compile: g77 -c -ffree-form ! Usage: N/A ! Arguments: N/A ! Requires: none ! Reference: sparspak.f90 from the web -- code of George & Liu 1981, see below !@Book{Geor:Liu:1981, !Author={George, A. and J.W.-H. Liu}, !Title={Computer solution of large sparse positive definite systems}, !Year=1981, !Address={Englewood Cliffs, New Jersey}, !Publisher={Prentice-Hall, Inc.}, !Pages=324 !} subroutine gsfct (n,ixlnz,lnz,xnzsub,nzsub,diag,link,first,temp,eps,war) ! !******************************************************************************* ! !! GS_FACTOR performs the symmetric factorization for a general sparse ! system, stored in the compressed subscript data format. ! ! Code taken from sparspak.f90 on the web (CJP 9/4/01) and amended to use ! Lockwood & Schervish generalized Cholesky factor (Lockwood et al. 2001) ! ! This is the double precision version ! ! Reference: ! ! Alan George and J W Liu, ! Computer Solution of Large Sparse Positive Definite Systems, ! Prentice Hall, 1981. ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! ixlnz-index vector for lnz. ixlnz(i) points to the ! start of nonzeros in column i of factor l. ! (xnzsub,nzsub)-the compressed subscript data ! structure for factor l. ! ! updated parameters- ! lnz-on input,contains nonzeros of a,and on ! return,the nonzeros of l. ! diag-the diagonal of l overwrites that of a. ! ! working parameters- ! link-at step j,the list in ! link(j),link(link(j)),. ! consists of those columns that will modify ! the column l(*,j). ! first-temporary vector to point to the first ! nonzero in each column that will be used ! next for modification. ! integer n ! CJP changed reals to double precision 9/26/01 double precision diag(n) double precision diagj integer first(n) integer i integer ii integer istop integer istrt integer isub integer ixlnz(n+1) integer j integer k integer kfirst integer link(n) double precision ljk double precision lnz(*) integer newk integer nzsub(*) double precision temp(n) integer xnzsub(*) double precision eps integer war integer writeout integer numcol ! ! Initialize working vectors. ! ! CJP changed these to for loops instead of single statements do 100 i= 1, n link(i) = 0 temp(i) = 0.0 100 continue ! ! Compute column l(*,j) for j = 1,....,N. ! numcol=0 do j = 1, n ! ! For each column l(*,k) that affects l(*,j). ! diagj = 0.0 newk = link(j) writeout=0 20 continue k = newk if ( k == 0 ) then goto 40 end if newk = link(k) ! ! Outer product modification of L(*,J) by L(*,K) starting at FIRST(K) of L(*,K). ! kfirst = first(k) ljk = lnz(kfirst) diagj = diagj + ljk*ljk if(writeout==1) then write(*,*) diagj, ljk end if istrt = kfirst+1 istop = ixlnz(k+1)-1 if ( istop < istrt ) then goto 20 end if ! ! Before modification, update vectors first, ! and link for future modification steps. ! first(k) = istrt i = xnzsub(k) + ( kfirst - ixlnz(k) ) + 1 isub = nzsub(i) link(k) = link(isub) link(isub) = k ! ! The actual mod is saved in vector temp. ! do ii = istrt, istop isub = nzsub(i) temp(isub) = temp(isub) + lnz(ii) * ljk if(writeout==1) then write( *,* ) temp(isub), lnz(ii), ljk end if i = i + 1 end do goto 20 ! ! Apply the modifications accumulated in temp to column L(*,J). ! 40 continue diagj = diag(j) - diagj if ( diagj <= eps ) then ! this commented out code will indicate where the problem with ! positive definiteness occurs ! if( war == 1) then ! war=j ! write ( *, * ) ' ' ! write ( *, * ) 'GSFCT - mat not pos def in column: ', j,diagj ! write ( *, * ) ' Zero or negative diagonal entry!' ! write ( *, * ) ' DIAG(J) = ', diagj ! write ( *, * ) ' for diagonal J = ', j ! write ( *, * ) ' ' ! write ( *, * ) ' The matrix is not positive definite.' ! end if ! ! this block of code added by CP 9/4/01 for generalized Cholesky numcol=numcol+1 diagj=0; diag(j)=0; istrt = ixlnz(j) istop = ixlnz(j+1) - 1 if ( istop >= istrt ) then first(j)=istrt i=xnzsub(j) isub=nzsub(i) link(j)=link(isub) link(isub)=j do ii = istrt, istop isub=nzsub(i) lnz(ii)=0 temp(isub)=0.0 i=i+1 end do end if goto 60 ! through here added by CP end if diagj = sqrt ( diagj ) diag(j) = diagj istrt = ixlnz(j) istop = ixlnz(j+1) - 1 if ( istop >= istrt ) then first(j) = istrt i = xnzsub(j) isub = nzsub(i) link(j) = link(isub) link(isub) = j do ii = istrt, istop isub = nzsub(i) lnz(ii) = ( lnz(ii) - temp(isub) ) / diagj if(writeout==1) then write(*,*) lnz(ii), temp(isub), diagj end if temp(isub) = 0.0 i = i+1 end do end if 60 end do ! if(numcol>0) then ! write(*,*) 'zeroed out columns:',numcol ! end if return end