C ****************************************************************** C ****************************************************************** subroutine moresor(n,g,bnnz,blin,bcol,bval,bdiag,delta,sigma1, +sigma2,eps,maxit,l,pd,p,chcnt,memfail,msinfo) implicit none C SCALAR ARGUMENTS logical memfail,pd integer bnnz,chcnt,maxit,msinfo,n double precision delta,eps,l,sigma1,sigma2 C ARRAY ARGUMENTS integer bcol(bnnz),bdiag(n),blin(bnnz) double precision bval(bnnz),g(n),p(n) C Solves the problem C C minimize psi(w) = 1/2 w^TBw + g^Tw C subject to ||w|| <= delta C C Using the method described in "Computing a trust region step", C by More and Sorensen. C msinfo: C C 0: both g and B are null C 1: first convergence criterion is satisfied C 2: second convergence criterion is satisfied C 3: third convergence criterion is satisfied C 5: maximum allowed number of iterations is achieved include "dim.par" include "machconst.inc" include "outtyp.inc" C LOCAL SCALARS integer col,dum,i,idx,iter,lin,lssinfo double precision b1n,d,delta2,geucn,ll,llant,ls,lsant,lu,luant, + peucn,peucn2,ptz,rpeucn2,rzeucn2,sgn,tau,teucn2,tmp,ueucn2 C LOCAL ARRAYS integer wi(nmax) double precision t(nmax),wd(nmax),z(nmax) delta2 = delta**2 msinfo = 0 memfail = .false. call lssini(.false.,.false.,.true.) C step 1: initialize ls (lower bound on l) with max{-bii}, where bii C are the elements of the diagonal of B ls = -bval(bdiag(1)) do i = 2,n lin = bdiag(i) ls = max( ls, -bval(lin) ) end do C Calculate ||B||1, B sparse do i = 1,n wd(i) = 0.0d0 end do do i = 1,bnnz lin = blin(i) col = bcol(i) if ( ( lin .le. n ) .and. ( col .le. n ) ) then wd(col) = wd(col) + abs( bval(i) ) if ( lin .ne. col ) then wd(lin) = wd(lin) + abs( bval(i) ) end if end if end do b1n = wd(1) do i = 2,n b1n = max( b1n, wd(i) ) end do C step 2: initialize ll (lower bound on l) with C max{0, ls, ||g||/delta - ||B||1}, where ||B||1 is the C 1-norm of the matrix B geucn = 0.0d0 do i = 1,n geucn = geucn + g(i)**2 end do geucn = sqrt( geucn ) ll = (geucn / delta) - b1n ll = max( 0.0d0, ll ) ll = max( ls, ll ) C step 3: initialize lu (upper bound on l) with ||g||/delta + ||B||1 lu = (geucn / delta) + b1n C If the matrix is null, there is nothing to be done if ( ( abs( ll ) .le. macheps23 ) .and. + ( abs( lu ) .le. macheps23 ) .and. + ( abs( ls ) .le. macheps23 ) ) then msinfo = 0 go to 21 end if C step 4: initialize iteration counter iter = 1 call lssana(n,bnnz,blin,bcol,bval,bdiag,lssinfo) if ( lssinfo .eq. 6 ) then ! INSUFFICIENT SPACE TO STORE THE LINEAR SYSTEM memfail = .true. return else if ( lssinfo .eq. 7 ) then ! INSUFFICIENT DOUBLE PRECISION WORKING SPACE memfail = .true. return else if ( lssinfo .eq. 8 ) then ! INSUFFICIENT INTEGER WORKING SPACE memfail = .true. return end if C step 5: safeguard of l (ensures that l is bigger than ll) 5 continue l = max( l, ll ) C step 6: safeguard of l (ensures that l is smaller than lu) l = min( l, lu ) C step 7: safeguard of l if ( ( l .le. ls + macheps23 * max( abs( ls ), 1.0d0 ) ) + .and. ( iter .ne. 1 ) ) then l = max( 1.0d-3 * lu, sqrt( ll*lu ) ) end if C step 8: try to use the Cholesky decomposition: (B +lI) = R^TR. C If the decomposition is successfull, R is stored in the C upper triangular portion of B (including the diagonal) and C pd is set to true. C If the decomposition fails, d and idx are set as explained C before, pd is set to false, and the Euclidian-norm of u is C calculated (see explanation of variable ueucn2) do i = 1,n wd(i) = l end do call lssfac(n,bnnz,blin,bcol,bval,bdiag,wd,idx,d,dum,lssinfo) chcnt = chcnt + 1 if ( lssinfo .eq. 6 ) then ! INSUFFICIENT SPACE TO STORE THE LINEAR SYSTEM memfail = .true. return else if ( lssinfo .eq. 7 ) then ! INSUFFICIENT DOUBLE PRECISION WORKING SPACE memfail = .true. return else if ( lssinfo .eq. 8 ) then ! INSUFFICIENT INTEGER WORKING SPACE memfail = .true. return end if if ( ( lssinfo .eq. 1 ) .or. ( lssinfo .eq. 2 ) ) then pd = .false. else ! lssinfo .eq. 0 pd = .true. end if C In this case (B + lI) is not positive definite, and d and idx are C calculated. Because p cannot be calculated (it is not possible to C solve the system using Cholesky factorization), the values of l, C ll and ls are updated, the iteration counter is increased and a C new iteration is started if ( .not. pd ) then C Print information (current iteration) if ( iprintinn .ge. 5 ) then C write(*, 1000) iter C write(*, 1010) ls,ll,lu,l C write(*, 1070) write(10,1000) iter write(10,1010) ls,ll,lu,l write(10,1070) end if llant = ll luant = lu lsant = ls if ( lu-l .le. macheps23 * max( abs( lu ),1.0d0 ) ) then lu = lu + macheps23 * max( abs( lu ),1.0d0 ) l = lu else call scalcu(n,bnnz,blin,bcol,bval,bdiag,l,idx,p,ueucn2,wd, + memfail) chcnt = chcnt + 1 if ( memfail ) then go to 22 end if ll = max( l, ll ) ls = max( l + (d / ueucn2), ls ) ll = max( ll, ls ) l = ls end if iter = iter + 1 C Test whether the number of iterations is exhausted if ( iter .gt. maxit ) then msinfo = 5 if ( iprintinn .ge. 5 ) then C write(*, 1090) write(10,1090) end if go to 22 end if go to 5 end if C step 9: solve R^TRp = -g for p and calculate the squared C Euclidian-norm of p do i = 1,n p(i) = -g(i) end do call lsssol(n,p) C Euclidian-norm of Rp = p^T R^TRp = p^T (-g) = - p^Tg rpeucn2 = 0.0d0 do i = 1,n rpeucn2 = rpeucn2 - p(i) * g(i) end do peucn2 = 0.0d0 do i = 1,n peucn2 = peucn2 + p(i)**2 end do peucn = sqrt( peucn2 ) C step 10: calculate z and tau, where tau * z is the approximation C of the eigenvector associated with the smallest C eigenvalue of B if ( peucn .lt. delta ) then C Calculate z call scalcz(n,wi,wd,z) C Calculate z Euclidian-norm tmp = 0.0d0 do i = 1,n tmp = tmp + z(i)**2 end do tmp = sqrt( tmp ) C Divide z by its norm do i = 1,n z(i) = z(i) / tmp end do C Calculate the squared Euclidian-norm of the product Rz. C Note that z^T R^T Rz = z^T (B + lI) z do i = 1,n wd(i) = 0.0d0 end do do i = 1,bnnz lin = blin(i) col = bcol(i) if ( lin .eq. col ) then wd(lin) = wd(lin) + (bval(i) + l) * z(col) else wd(lin) = wd(lin) + bval(i) * z(col) wd(col) = wd(col) + bval(i) * z(lin) end if end do rzeucn2 = 0.0d0 do i = 1,n rzeucn2 = rzeucn2 + z(i) * wd(i) end do C Calculate tau ptz = 0.0d0 do i = 1,n ptz = ptz + p(i) * z(i) end do if ( ptz .lt. 0.0d0 ) then sgn = -1.0d0 else sgn = 1.0d0 end if tmp = delta2 - peucn2 tau = (ptz**2) + tmp tau = tmp / (ptz + sgn * sqrt( tau )) end if C Print informations (current iteration) if ( iprintinn .ge. 5 ) then C write(*, 1000) iter C write(*, 1010) ls,ll,lu,l C write(*, 1020) peucn write(10,1000) iter write(10,1010) ls,ll,lu,l write(10,1020) peucn if ( peucn .lt. delta ) then C write(*, 1030) sqrt( peucn2 + 2 * tau * ptz + tau ** 2 ) write(10,1030) sqrt( peucn2 + 2 * tau * ptz + tau ** 2 ) else C write(*, 1030) peucn write(10,1030) peucn end if C write(*, 1050) delta write(10,1050) delta end if C steps 11 and 12: update ll, lu and ls llant = ll luant = lu lsant = ls if ( peucn .lt. delta ) then lu = min( l, lu ) ls = max( l - rzeucn2, ls ) else ll = max( l, ll ) end if C step 13: update ls when B + lI is not positive definite. C This was done right after the Cholesky decomposition C failure C step 14: update ll ll = max( ll, ls ) C step 15: convergence test if ( ( abs( l ) .le. eps ) .and. ( peucn .le. delta ) ) then msinfo = 1 go to 21 end if C step 16: second convergence test if ( abs( delta - peucn ) .le. sigma1*delta ) then msinfo = 2 end if C step 17: convergence test for the hard case tmp = rpeucn2 + l * delta2 tmp = max( sigma2, tmp ) tmp = tmp * sigma1 * (2 - sigma1) if ( ( peucn .lt. delta ) .and. ( (rzeucn2*(tau**2) .le. tmp ) + .or. (lu - ll .le. eps) ) ) then msinfo = msinfo + 3 go to 20 end if if ( msinfo .eq. 2 ) then go to 21 end if C step 21: Calculate l to be used in the next iteration if ( ( abs( geucn ) .gt. eps ) .or. + ( ( l .le. ll + macheps23 * max( abs( ll ), 1.0d0 ) ) + .and. ( ls .le. ll ) ) ) then C Solve R^T t = p for t and calculate t squared Euclidian-norm do i = 1,n t(i) = p(i) end do call lsssoltr('T',n,t) teucn2 = 0.0d0 do i = 1,n teucn2 = teucn2 + t(i)**2 end do C Update l using Newton's method update l = l + (peucn2 / teucn2) * ((peucn - delta) / delta) else l = ls end if C step 22: update iteration counter iter = iter + 1 C Test whether the number of iterations is exhausted if ( iter .gt. maxit ) then msinfo = 5 if ( iprintinn .ge. 5 ) then C write(*, 1090) write(10,1090) end if go to 22 end if C step 23: start a new iteration if ( ( abs( llant-ll ) .le. + macheps23 * max( abs( ll ), 1.0d0 ) ) .and. + ( abs( luant-lu ) .le. + macheps23 * max( abs( lu ), 1.0d0 ) ) .and. + ( abs( lsant-ls ) .le. + macheps23 * max( abs( ls ), 1.0d0 ) ) ) then ll = ll + macheps23 * max( abs( ll ), 1.0d0 ) lu = lu - macheps23 * max( abs( lu ), 1.0d0 ) end if go to 5 C steps 18, 19 and 20: C The solution is given by p in 3 cases: C - if the first convergence criterion is satisfied; C - if only the second convergence criterion is satisfied; C - if both the second and the third convergence criteria are C satisfied, but the squared Euclidian-norm of R(tau*z) is C strictly bigger than l*(delta2 - peucn2) C The solution is given by p + tau*z when: C - just the third convergence criterion is satisfied; C - both the second and the third convergence criteria are C satisfied, but the squared Euclidian-norm of R(tau*z) is smaller C or equal to l*(delta2 - peucn2) 20 tmp = (rzeucn2 * (tau**2)) - l * (delta2 - peucn2) if ( ( msinfo .eq. 3 ) .or. + ( ( msinfo .eq. 5 ) .and. ( tmp .le. 0.0d0 ) ) ) then peucn2 = 0.0d0 do i = 1,n p(i) = p(i) + tau * z(i) peucn2 = peucn2 + p(i)**2 end do peucn = sqrt( peucn2 ) msinfo = 3 else msinfo = 2 end if C Print informations 21 if ( iprintinn .ge. 5 ) then C write(*, 1060) iter C write(*, 1010) ls,ll,lu,l C write(*, 1080) C write(*, 1040) peucn write(10,1060) iter write(10,1010) ls,ll,lu,l write(10,1080) write(10,1040) peucn end if 22 continue C Non-executable statements 1000 format(/,10X,'More-Sorensen iteration: ',I7) 1010 format( 10X,'ls = ',1P,D11.4, + 10X,'ll = ',1P,D11.4, + 10X,'lu = ',1P,D11.4, + 10X,'l = ',1P,D11.4) 1020 format( 10X,'Euclidian-norm of p: ',1P,D7.1) 1030 format( 10X,'Euclidian-norm of p + z: ',1P,D7.1) 1040 format( 10X,'Euclidian-norm of step: ',1P,D7.1) 1050 format( 10X,'delta: ',1P,D7.1) 1060 format( 10X,'Number of iterations: ',I7) 1070 format( 10X,'Matrix is not positive definite!') 1080 format( 10X,'Flag of More-Sorensen: ', + 'convergence criterion satisfied.') 1090 format( 10X,'Flag of More-Sorensen: ', + 'maximum number of iterations achieved.') end C ****************************************************************** C ****************************************************************** subroutine scalcu(n,annz,alin,acol,aval,adiag,l,idx,u,ueucn2,wd, +memfail) implicit none C SCALAR ARGUMENTS logical memfail integer annz,idx,n double precision l,ueucn2 C ARRAY ARGUMENTS integer acol(annz),adiag(n),alin(annz) double precision aval(annz),u(idx),wd(idx-1) C Solves a sparse linear system of the form (A + lI + ekek^Td)u = 0, C where A + lI is a definite positive matrix in R^{k x k}. C u is a vector of k positions with u(k) = 1. C LOCAL SCALARS integer col,i,lin,lssinfo if ( idx .eq. 1 ) then u(idx) = 1.0d0 ueucn2 = 1.0d0 return end if do i = 1,idx u(i) = 0.0d0 end do C Permute columns and rows of A call lsspermind(annz,alin,acol) C Eliminate columns that have index greater than idx and define u C as idx-th column of A + lI do i = 1,annz col = acol(i) lin = alin(i) if ( ( col .eq. idx ) .and. ( lin .lt. idx ) ) then u(lin) = u(lin) - aval(i) else if ( ( lin .eq. idx ) .and. ( col .lt. idx ) ) then u(col) = u(col) - aval(i) end if if ( ( col .ge. idx ) .or. ( lin .ge. idx ) ) then acol(i) = col + n alin(i) = lin + n end if end do C Solve system (A + lI)x = u call lssini(.false.,.true.,.false.) call lsspermvec(n,adiag) do i = 1,idx-1 wd(i) = l end do call lssafsol(idx-1,annz,alin,acol,aval,adiag,wd,u,lssinfo) if ( lssinfo .eq. 6 ) then ! INSUFFICIENT SPACE TO STORE THE LINEAR SYSTEM memfail = .true. return else if ( lssinfo .eq. 7 ) then ! INSUFFICIENT DOUBLE PRECISION WORKING SPACE memfail = .true. return else if ( lssinfo .eq. 8 ) then ! INSUFFICIENT INTEGER WORKING SPACE memfail = .true. return end if call lssunpermvec(n,adiag) call lssini(.false.,.false.,.true.) C Undo modifications in acol and alin do i = 1,annz col = acol(i) if ( col .gt. n ) then acol(i) = acol(i) - n alin(i) = alin(i) - n end if end do call lssunpermind(n,annz,alin,acol) u(idx) = 1.0d0 ueucn2 = 0.0d0 do i = 1,idx ueucn2 = ueucn2 + u(i)**2 end do end C ****************************************************************** C ****************************************************************** subroutine scalcz(n,rowind,rowval,z) implicit none C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS integer rowind(n) double precision rowval(n),z(n) C Sparse implementation of the technique presented by Cline, Moler, C Stewart and Wilkinson to estimate the condition number of a upper C triangular matrix. include "machconst.inc" C LOCAL SCALARS integer col,i,k,rownnz double precision acs,acz,d,ek,rki,s,sm,w,wk,wkm double precision lssgetd C EXTERNAL FUNCTIONS external lssgetd call lsssetrow(n) do i = 1,n z(i) = 0.0d0 end do acz = 0.0d0 acs = 1.0d0 ek = 1.0d0 do k = 1,n d = lssgetd(k) if ( abs( acs*z(k) ) .gt. macheps23 ) then ek = dsign(ek,-acs*z(k)) end if if ( abs( ek - acs*z(k) ) .gt. abs( d ) ) then s = abs( d ) / abs( ek - acs*z(k) ) acs = acs * s ek = s * ek end if wk = ek - acs * z(k) wkm = -ek - acs * z(k) s = abs( wk ) sm = abs( wkm ) if ( d .eq. 0.0d0 ) then wk = 1.0d0 wkm = 1.0d0 else wk = wk / d wkm = wkm / d end if if ( k .eq. n ) then go to 10 end if call lssgetrow(n,k,rownnz,rowind,rowval) sm = sm + acs * acz do i = 1,rownnz col = rowind(i) if ( col .gt. k ) then rki = rowval(i) sm = sm - abs( acs * z(col) ) sm = sm + abs( acs * z(col) + wkm * rki ) acz = acz - abs( z(col) ) acz = acz + abs( z(col) + (wk/acs) * rki ) z(col) = z(col) + (wk/acs) * rki end if end do s = s + acs * acz if ( s .lt. + sm - macheps23 * max( abs( sm ), 1.0d0 ) ) then w = wkm - wk wk = wkm do i = 1,rownnz col = rowind(i) if ( col .gt. k ) then rki = rowval(i) acz = acz - abs( z(col) ) acz = acz + abs( z(col) + (w/acs) * rki ) z(col) = z(col) + (w/acs) * rki end if end do end if acz = acz - abs( z(k+1) ) 10 z(k) = wk/acs end do C Divide z by its 1-norm to avoid overflow s = 0.0d0 do i = 1,n s = s + abs( z(i) ) end do do i = 1,n z(i) = z(i) / s end do C Solve Rz = y call lsssoltr(' ',n,z) end C ****************************************************************** C ****************************************************************** C moresor: C C Method to minimize sparse quadratic functions subjected to C ||w|| <= delta C C minimize psi(w) = 1/2 w^TBw + g^Tw C subject to ||w|| <= delta C C Method described in "Computing a trust region step", by More and C Sorensen. C C The main ideia of this method is to find a positive scalar \mslamb C that is a zero of the function C C phi(\mslamb) = 1/||p|| - 1/delta, C C where p is the solution of the linear system C C (B + \mslamb I)p = -g. C C Note that symmetric matrix B, vector g and positive real number C delta are presented in the minimization problem above. I is the C identity matrix. C C The method used to find the zero of that function is basically the C Newton method to find roots. C C On Entry C C n integer C dimension C C g double precision g(n) C vector used to define the quadratic function C C bnnz integer C number of nonzero elements of B C C blin integer blin(bnnz) C row indices of nonzero elements of B C C bcol integer bcol(bnnz) C column indices of nonzero elements of B C C bval double precision bval(bnnz) C nonzero elements of B, that is, C B(blin(i),bcol(i)) = bval(i). Since B is symmetric, just C one element B(i,j) or B(j,i) must appear in its sparse C representation. If more than one pair corresponding to C the same position of B appears in the sparse C representation, the multiple entries will be summed. C If any blin(i) or bcol(i) is out of range, the entry will C be ignored C C bdiag integer bdiag(n) C indices of diagonal elements of B in blin, bcol and bval C C delta double precision C trust-region radius C C sigma1 double precision C allowed error for convergence criteria 1 and 2 C C sigma2 double precision C allowed error for convergence criteria 3 (hard case) C C eps double precision C allowed error C C maxit integer C maximum number of allowed iterations C C l double precision C initial value for mslamb C C On Return C C l double precision C value that gives p as a solution to the minimization C problem, because p is also solution to C (B + l I)p = -g C C pd logical C set to true if the last Cholesky decomposition is C successfull C C p double precision p(n) C solution to problem C minimize psi(w) C subjected to ||w|| <= delta C C chcnt integer C number of Cholesky decompositions C C memfail logical C true iff linear solver failed because of lack of memory C C msinfo integer C stores which convergence criteria was satisfied: C C 0 = both g and B are null; C C 1 = first convergence criterion is satisfied; C C 2 = second convergence criterion is satisfied; C C 3 = third convergence criterion is satisfied C C 5 = maximum allowed number of iterations is achieved. C ****************************************************************** C ****************************************************************** C scalcu: C C Solve a sparse linear system of the form (A + lI + ekek^Td)u = 0, C where A + lI is a definite positive matrix in R^{k x k}. C u is a vector of k positions with u(k) = 1. C C On Entry C C n integer C dimension of A C C annz integer C number of nonzero elements of A C C alin integer alin(annz) C row indices of nonzero elements of A C C acol integer acol(annz) C column indices of nonzero elements of A C C aval integer aval(annz) C nonzero elements of A, that is, C A(alin(i),acol(i)) = aval(i). Since A is symmetric, just C one element A(i,j) or A(j,i) must appear in its sparse C representation. If more than one pair corresponding to C the same position of A appears in the sparse C representation, the multiple entries will be summed. C If any alin(i) or acol(i) is out of range, the entry will C be ignored C C adiag integer adiag(n) C indices of diagonal elements of A in alin, acol and aval C C l double precision C used to compute A + lI C C idx integer C index k C C On Return C C u double precision u(n) C system solution C C ueucn2 double precision C u squared Euclidian-norm C C memfail logical C true iff linear solver failed because of lack of memory C ****************************************************************** C ****************************************************************** C scalcz: C C Sparse implementation of the technique presented by Cline, Moler, C Stewart and Wilkinson to estimate the condition number of a matrix. C This technique is used by the More-Sorensen method to calculate an C approximation to the eigenvector associated to the smallest C eigenvalue of a matrix B (\mslamb_1). C In this technique, when \mslamb approaches -\mslamb_1, \|Rz\| C approaches 0. This insures that z is an approximation to the C wanted eigenvector. C Basically, it solves R^Ty = e, choosing e(k) as 1 or -1 (whatever C gives maximum local growth of y). Then, it solves Rz = y. C Note that R is a sparse matrix given by P^T D^0.5 L^T P (obtained C applying subroutine MA27BD). C C On Entry C C n integer C dimension of R C C rowind integer rowind(n) C rowval double precision rowval(n) C working arrays C C On Return C C z double precision z(n) C approximation of the eigenvector of B associated to C \mslamb_1