! File: gsfct.real.f ! Date: 12/31/01 (downloaded fall 2001) ! Author: C. Paciorek - taken off the web ! Purpose: Cholesky factorization of pos def matrix stored in sparse format ! single precision version ! Compile: g77 -c -ffree-form ! Usage: ! Arguments: ! Requires: ! Reference: sparspak.f90 from the web -- George & Liu 1981, see below subroutine gsfct (n,ixlnz,lnz,xnzsub,nzsub,diag,link,first,temp,epsilon) ! !******************************************************************************* ! !! 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 ! Schervish' generalized Choleski factor (Lockwood et al.) ! ! 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 real diag(n) real 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) real ljk real lnz(*) integer newk integer nzsub(*) real temp(n) integer xnzsub(*) real epsilon ! ! 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. ! do j = 1, n ! ! For each column l(*,k) that affects l(*,j). ! diagj = 0.0 newk = link(j) 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 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 i = i + 1 end do goto 20 ! ! Apply the modifications accumulated in temp to column L(*,J). ! 40 continue diagj = diag(j) - diagj ! JR used .000001 if ( diagj <= epsilon ) then ! write ( *, * ) ' ' ! write ( *, * ) 'GSFCT-matrix not pos def in col: ', 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.' ! diagj=0; ! this block of code added by CP 9/4/01 for gen chol 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) ! write (*,*) 'GSF ', lnz(ii) , temp(isub) , diagj lnz(ii) = ( lnz(ii) - temp(isub) ) / diagj ! write (*,*) 'sf ', lnz(ii) , temp(isub) , diagj temp(isub) = 0.0 i = i+1 end do end if 60 end do return end