!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! PROGRAM VPM (Variable Projection Method) ! ! AUTHORS: Valeria Ruggiero (Department of Mathematics - ! University of Ferrara) ! Luca Zanni (Department of Mathematics - University ! of Modena and Reggio Emilia) ! ! 23/09/99 ! ! ! This program computes by VPM a solution of a linearly ! constrained convex quadratic programming (QP) problem of ! the form ! ! minimize f(x)=1/2 xt A x + qt x ! subject to C1 x >=b1 ! C2 x =b2 ! ! where x is the n-dimensional unknown vector; ! m is the number of all the contraints; A is an n*n matrix; ! C1 is an (m-meq)*n matrix; C2 is an neq*n matrix of full ! row-rank and ! _ _ _ _ ! | C1 | | b1 | ! C = | | b=| | ! |_ C2 _| |_b2_| ! ! When the QP problem has inequality constraints, the sequence of ! subproblems generated by the VPM are solved by a projected SOR ! method (Cryer, SIAM J. Control, 1971). ! In presence of equality constraints only, VPM uses a ! Preconditioned Conjugate Gradient Method (PCG) as inner solver. ! This last algorithm uses the Arithmetic Mean Preconditioner. ! (Galligani - Ruggiero, International Journal of Computer ! Mathematics, 1992). ! In input, the user has to introduce the data described below ! (see "Description of the input variable" and the file inputvpm ! that is an example of data input). ! ! The program uses double precision floating point arithmetic: in ! the description, the word "real" in the description of input ! data stands for double precision real. ! ! For some basic linear algebra operation, the program uses the ! BLAS1 routines. ! ! We point out the way the sparse matrices A and B are stored: each row contains ! only the nonzero elements of A (or B), then a vector ka (or kb) is used to keep the ! number of non-zero elements of each row and a matrix ia (or ib) is arranged so that ! ia(i,*) contains the column index of the elements in the i-th row of A. ! EXAMPLE: ! _ _ _ _ _ _ _ _ ! | 3 0 2 | | 2 | | 1 3 | | 3 2 | ! A= | 0 4 1 | ==> ka= | 2 | ia=| 2 3 | a=| 4 1 | ! |_ 2 1 4_| |_ 3_| |_1 2 3_| |_2 1 4_| ! !--------------- DESCRIPTION OF THE INPUT VARIABLES ---------------------------------- ! ! name = name for the file containing the matrix A (character) ! namev = name for the file containing the constraints (character) ! omega = relaxation parameter for the projected SOR method or the ! Arithmetic Mean Preconditioner (real) ! eps = tolerance for the stopping criterion of VPM ! (outer tolerance) (real) ! epstart = maximum tolerance for the stopping criterion of projected ! SOR method (inner solver); the projected SOR method uses ! an inner progressive termination rule, depending on the ! quality of the previous iterate; if the PCG is used as ! inner solver, epstart is ignored; (real) ! epssor = minimum tolerance for the the stopping criterion of ! projected SOR method (inner solver); if the PCG is used ! as inner solver epssor is used in the inner stopping ! rule; (real) ! indrule = if indrule=1, the updating of the projection parameter ! following the rule ! (dt A d)/(dt A inv(D) A d) d=descent direction ! otherwise we use the rule ! (dt D d)/(dt A d) (integer) ! ! The unformatted files "name" and "namev" contain the data of the QP ! problem as described below. ! ! FILE "name" contains the following records: ! - n, inc, nrango = order of A, number of nonzero elements of the ! upper triangular part of A, rank of A;(integer) ! - A = matrix A; A is stored in compressed form, as a set of ! consecutive records which contain the nonzero elements of ! the upper triangular part of A, in a row wise order. ! Each record has the informations of four nonzero elements ! of A in the following format: ! (row index i, column index j, element (i,j) of the matrix) ! - q = vector of n element (real) ! - fmin= minimum value of the objective function f(x) (real) ! if fmin is unknown, set fmin=1d-25; ! - sparsa= sparsity of A (percentage) (real) ! FILE "namev" contains the following records: ! - m,meq,nav,inc = number of constraints, number of equality ! constraints, ! number of active inequality constraints in a ! solution ! number of nonzero elements of the matrix of the ! constraints; (integer) ! - C = matrix C; C is stored in compressed form, with the same ! format of matrix A; all the nonzero elements of B are ! given; each record has the informations of four nonzero ! elements of C in the following format: ! (row index i, column index j, element (i,j) of the matrix) ! - xe = vector containing a solution of the QP problem (real) ! if xe is unknown, set xe(i)=1d-25 for any i; ! - b = vector of m elements (real) ! - sparb = sparsity of C (percentage) (real) ! ! The files name and namev have the same format of those generated by ! program TEST98. !------------------------------------------------------------------------ !---------- OUTPUT: DESCRIPTION ------------------------------------------ ! ! As output, the program gives a formatted file "vpm.dat" ! containing a set of input informations (about the input data ! and the features of the problem) and the output data, ! computed by the routine VPMSOR or VPMPCG and by the program. ! In particular, the routine VPMSOR or VPMPCG gives the ! following results: ! ! info = if info=1, VPM does not converge in maxit ! outer iterations; if info=2, the inner solver has ! performed the maximum number of iterations at least ! one time; ! iter = number of VPM iterations ! icont = total number of inner iterations ! nacc = number of outer iterations using the correction rule ! X = vector of the computed solution; ! YZ = vector of Lagrange multipliers associated to the computed ! solution; ! ! Furthermore, the program gives the following informations: ! ! fmin1 = value of f(x) computed at the obtained solution ! fnorm = relative error of the minimum value of f(x) ! (if a correct value of fmin is given in input) ! xnorm = relative errore of the solution (2-norm) ! (if the exact vector solution of the problem ! is given in input) ! nac = number of active constraints active in the computed ! solution ! t2 = elapsed time in seconds ! ! !******* REFERENCES **************************************************** ! ! 1. V. Ruggiero, L. Zanni - On the efficiency of splitting and ! projection methods for large strictly convex quadratic ! programs - "Nonlinear Optimization and Related Topics", ! (G. Di Pillo, F. Giannessi eds.), ! Kluwer Academic Publishers (1999). ! 2. V. Ruggiero, L. Zanni - A modified projection algorithm for ! large strictly convex quadratic programs - Journal of ! Optimization Theory and Applications, Vol.104, n. 2 (2000). ! 3. V. Ruggiero L. Zanni - An Overview on Projection-Type ! Methods for Convex Large-Scale Quadratic Programs - "Equilibrium ! Problems and Variational Models", (A. Maugeri, F. Giannessi, ! P. Pardalos eds.), Kluwer Academic Press (2000). ! 4. V. Ruggiero L. Zanni -Variable Projection Methods for Large ! Convex Quadratic Programs - "Recent Trend in Numerical Analysis", ! (L. Brugnano, D. Trigiante eds.), Advances in Computation Theory ! and Practice, Nova Science Book and journals, Huntington, NJ (2000). ! !***************************************************************************************** PROGRAM VPM IMPLICIT DOUBLE PRECISION(a-h,o-z) PARAMETER(in=5001,im=5001,ibmax=300,icmax=100, 1 maxit=5000,maxit1=3000, 2 inm=icmax*im,inq=ibmax*in) ! ! in = maximum number of rows of matrix A ! im = maximum number of rows of matrix C ! ibmax = maximum number of nonzero elements belonging to a row ! of A ! icmax = maximum number of nonzero elements belonging to a row ! of C ! maxit = maximum number of outer iterations (VPM) ! maxit1= maximum number of inner iterations (projected SOR ! or PCG) ! ! real secnds ! ! The FUNCTION SECNDS(t1) for measuring the elapsed time from t1 ! is processor dependent: its substitution with a suitable one ! is left to the user. t1 is a real parameter. ! COMMON/DATA/ C(im,icmax),B(im),Q(in),A(in,ibmax),XE(in), 1IA(in,ibmax),IC(im,icmax),K1(in),K2(in),KC(im) COMMON/DIAGO/ X(0:in),YZ(im),S(im),S1(0:in), 1YZ0(im),HX(in),DD(IN) DATA C/inm*0.0/,A/inq*0.0/,IC/inm*0/,IA/inq*0/ DATA B/im*0.0/,Q/in*0.0/ ! CALL leggif(n,iban,m,meq,nva,ibanc,fmin,nrango, 1 sparsa,sparsc) ! READ(*,*) OMEGA IF (OMEGA.eq.1.0) then OMEGA=1.1 ENDIF READ(*,*) eps, epstart, epssor READ(*,*) indrule ! ! Compute the norm of A ! DO i=1,n X(i)=0d0 DO j=1,iban X(i)=X(i)+ABS(A(i,j)) ENDDO ENDDO RELE3=ABS(X(IDAMAX(n,X(1),1))) ! ! starting vector (this vector is, in general, not feasible) ! DO i=1,n X(i)=0.0 ENDDO ! ! choice of diagonal matrix ! do i=1,n DD(i)=1. enddo ! ! gamma= inverse of the projection parameter ! ! initial value for gamma gamma=0d0 DO i=1,n DO j=1,iban gamma=gamma+A(i,j)**2 ENDDO ENDDO gamma=sqrt(gamma) /DD(idamin(n,DD,1)) open(unit=1,file='vpm.dat',form='formatted') WRITE(1,*) WRITE(1,*) 'Variable projection Method' IF(m-meq.eq.0)then WRITE(1,*)'combined as PCG as inner solver' ELSE WRITE(1,*)'combined as projected SOR as inner solver' ENDIF WRITE(1,*) WRITE(1,*)'Diagonal matrix D=identity' WRITE(1,*) WRITE(1,*)'***INPUT PARAMETERS***' WRITE(1,*)'size of the problem =',n WRITE(1,*)'number of the linear constraints =', m WRITE(1,*)'number of the equality constraints=',meq WRITE(1,*)'rank of A=',nrango WRITE(1,*)'norm of A=',rele3 WRITE(1,*)'sparsity of A (percentage)=',sparsa WRITE(1,*)'sparsity of C (percentage)=',sparsc WRITE(1,*)'maximum number of nonzero elements per row' WRITE(1,*)'of A =',iban WRITE(1,*)'of C =',ibanc WRITE(1,*)'updating rule for the projection parameter=',indrule WRITE(1,*)'initial projection parameter=', 1/gamma WRITE(1,*)'outer tolerance =', eps IF(m-meq.eq.0)then WRITE(1,*)'inner tolerance (PCG) =',epssor ELSE WRITE(1,*)'minimum inner tolerance (epstart)=',epstart WRITE(1,*)'maximum inner tolerance (epssor) =',epssor WRITE(1,*)'inner tolerance: ' WRITE(1,*)'epsk=min(epstart,max(kkt**1.5,epssor))' WRITE(1,*)'where kkt= norm of the vector A x + q - Bt yz' ENDIF WRITE(1,*)'relaxation parameter =',omega WRITE(1,*) IF(m-meq.ne.0)then t1=secnds(real(0.)) CALL VPMSOR(n,m,meq,eps,maxit,maxit1,omega,info,iter, 1icont,iban,ibanc,kirmax,epssor,epstart,rele3,gamma, 2nacc,indrule) t2=secnds(real(t1)) ELSE t1=secnds(real(0.)) CALL VPMPCG(n,m,meq,eps,maxit,maxit1,omega,info,iter, 1icont,iban,ibanc,kirmax,epssor,epstart,rele3,gamma, 2nacc,indrule) t2=secnds(real(t1)) ENDIF IF(info.eq.1) then WRITE(1,*) 'VPM does not converge' STOP ENDIF WRITE(1,*)'info =',info WRITE(1,*)'Elapsed time in seconds =',t2 WRITE(1,*)'Number of VPM iterations =', iter WRITE(1,*)'Total number of inner iterations=', icont WRITE(1,*)'Maximum number of nonzero elements per row' WRITE(1,*)'in the matrix of linear complementarity ' WRITE(1,*)'problem (or in the inner linear system =',kirmax WRITE(1,*)'number of outer iterations using the ' WRITE(1,*)'correction rule = ',nacc ! ! Compute the objective function at the obtained numerical solution ! fmin1=0.0d0 DO i=1,n S1(i)=0.0d0 DO j=1,iban S1(i)=S1(i)+A(i,j)*X(IA(i,j)) ENDDO ENDDO DO i=1,n fmin1=fmin1+S1(i)*X(i) ENDDO fmin1=fmin1*0.5 DO i=1,n fmin1=fmin1+Q(i)*X(i) ENDDO ! ! Compute the relative error of the minimum value of f(x) ! WRITE(1,*) 'Computed minimum value of f(x) = ',fmin1 IF(fmin.ne.1d-25)THEN fnorm=abs(fmin-fmin1)/abs(fmin) WRITE(1,*) 'Exact minimum value of f(x) = ',fmin WRITE(1,*) 'Relative error = ',fnorm WRITE(1,*) ENDIF ! ! Compute the relative error of the solution in 2-norm ! DO i=1,n IF(xe(I).ne.1d-25)goto 2000 ENDDO GOTO 3000 2000 xnorm=dnrm2(n,xe,1) DO i=1,n XE(i)=XE(i)-X(i) ENDDO xnorm=DNRM2(n,XE,1)/xnorm WRITE(1,*)'relative error of the solution (2-norm) =',xnorm 3000 continue ! ! Compute the active constraints in the solution ! DO i=1,m-meq XE(i)=B(i) DO j=1,ibanc XE(i)=XE(i)-C(i,j)*X(IC(i,j)) ENDDO ENDDO nac=meq DO i=1,m-meq IF(ABS(XE(i)) .le. 5*eps) nac=nac+1 ENDDO WRITE(1,*)'number of active constraints = ',nac WRITE(1,*) WRITE(1,*)'********Solution***********' WRITE(1,'(1x,i4,7x,e14.7)')(i, x(i),i=1,n) WRITE(1,*)'********Vector of Lagrance multipliers***********' WRITE(1,'(1x,i4,7x,e14.7)')(i,yz(i),i=1,m) WRITE(1,*) '***********************' STOP END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Subroutine for reading the input data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE leggif(n,iban,mm,meq,nva,ibanc,fmin,nrango,sparsa, 1 sparsc) IMPLICIT DOUBLE PRECISION(a-h,o-z) PARAMETER(in=5001,im=5001,ibmax=300,icmax=100) COMMON/DATA/ C(im,icmax),B(im),Q(in),A(in,ibmax),AK(in), 1IA(in,ibmax),IC(im,icmax),K1(in),K2(in),KC(im) DIMENSION l(4),m(4),coeff(4) CHARACTER*15 name,namev ! ! input of file names ! READ(*,'(a15)') name,namev ! ! file of objective function ! OPEN(unit=8,file=name,form='unformatted',status='old') READ(8) n,inc,nrango IF(n.gt.in)THEN WRITE(*,*)'error: IN is small!!' STOP ENDIF DO i=1,n IA(i,1)=i K2(i)=0 K1(i)=1 ENDDO DO i=1,inc/4 READ(8) (l(j),m(j),coeff(j),j=1,4) DO j=1,4 IF(l(j).eq.m(j))THEN A(l(j),1)=coeff(j) ELSE K1(l(j))=K1(l(j))+1 IF(K1(l(j)).gt.ibmax)THEN WRITE(*,*)'error: IBMAX is small!!' STOP ENDIF A(l(j),K1(l(j)))=coeff(j) IA(l(j),k1(l(j)))=m(j) ENDIF ENDDO ENDDO ! irest=mod(inc,4) IF(irest. ne. 0)THEN READ(8)(l(j),m(j),coeff(j),j=1,irest) DO j=1,irest IF(l(j).eq.m(j))THEN A(l(j),1)=coeff(j) ELSE K1(l(j))=K1(l(j))+1 IF(K1(l(j)).gt.ibmax)THEN WRITE(*,*)'error: IBMAX is small!!' STOP ENDIF A(l(j),K1(l(j)))=coeff(j) IA(l(j),K1(l(j)))=m(j) ENDIF ENDDO ENDIF READ(8) (Q(i),i=1,n) READ(8) fmin READ(8) sparsa CLOSE(8) DO i=1,n DO j=2,K1(i) K2(IA(i,j))=K2(IA(i,j))+1 icol=K1(IA(i,j))+K2(ia(i,j)) IF(icol.gt.ibmax)THEN WRITE(*,*)'error: IBMAX is small!!' STOP ENDIF A(IA(i,j),icol)=A(i,j) IA(IA(i,j),icol)=i ENDDO ENDDO DO i=1,n AK(i)=K1(i)+K2(i) ENDDO iban=AK(idamax(n,AK,1)) ! ! file of data constraints ! OPEN(unit=8,file=namev,form='unformatted',status='old') READ(8) mm,meq,nva,inc IF(MM.gt.im)then WRITE(*,*)'error: IM is small!!' STOP ENDIF DO i=1,mm KC(i)=0 ENDDO DO i=1,inc/4 READ(8) (l(j),m(j),coeff(j),j=1,4) DO j=1,4 KC(l(j))=KC(l(j))+1 IF(KC(l(j)).gt.icmax)THEN WRITE(*,*)'error: ICMAX is small!!' STOP ENDIF C(l(j),KC(l(j)))=coeff(j) IC(l(j),KC(l(j)))=m(j) ENDDO ENDDO irest=mod(inc,4) IF(irest .ne. 0) THEN READ(8)(l(j),m(j),coeff(j),j=1,irest) DO j=1,irest KC(l(j))=KC(l(j))+1 IF(KC(l(j)).gt.icmax)THEN WRITE(*,*)'error: ICMAX is small!!' STOP ENDIF C(l(j),KC(l(j)))=coeff(j) IC(l(j),KC(l(j)))=m(j) ENDDO ENDIF DO i=1,mm B(i)=KC(i) ENDDO ibanc=B(idamax(mm,B,1)) READ(8)(AK(j),j=1,n) READ(8)(B(j),j=1,mm) READ(8) sparsc RETURN END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Subroutine of VPM combined with the projected SOR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE VPMSOR(n,m,meq,eps,maxit,maxit1,omega, 1info,iter,icont,iban,ibanc,kirmax,epssor,epstart, 2rele3,gamma,nacc,indrule) ! IMPLICIT DOUBLE PRECISION(a-h,o-z) PARAMETER (in=5001,im=5001,ibmax=300,icmax=100, 1 tol=1d-30,icc=700,iw=im+in) ! ! tol = the element of the complementarity matrix ! less than tol are considered as zero ! icc = maximum number of nonzero elements per row ! of the complementarity matrix ! DIMENSION CM(im,icc),ICM(im,icc),AINV(in),HXOLD(in), 1 SN(iw),TEMP(in),X0(in),IR(im),DIA(im),TTT(0:in) COMMON/DATA/ C(im,icmax),B(im),Q(in),A(in,ibmax),XE(in), 1IA(in,ibmax),IC(im,icmax),K1(in),K2(in),KC(im) COMMON/DIAGO/ X(0:in),YZ(im),S(im),S1(0:in), 1YZ0(im),HX(in),DD(in) DATA IR/im*1/ ! ! info=0 S1(0)=0.0 X(0)=0.0 TTT(0)=0.0 DO i=1,m YZ(i)=0.0 ENDDO DO i=1,n X0(i)=X(i) AINV(i)=1./DD(I) ENDDO mdq=m-meq epsk=epstart ! ! compute the complementarity matrix; this matrix is ! stored in compressed form in CM (ICM and KC) ! DO i=1,m DO j=1,n TEMP(j)=0.0 ENDDO sx=0.0 DO j=1,kc(i) TEMP(IC(i,j))=C(i,j)*AINV(IC(i,j)) sx=sx+TEMP(IC(i,j)) * C(i,j) ENDDO CM(i,1)=sx ICM(i,1)=i DO j=i+1,m sx=0. DO l=1,KC(j) sx=sx+C(j,l)*TEMP(IC(j,l)) ENDDO IF(ABS(sx) .gt. tol)THEN IR(i)=IR(i)+1 IF(IR(i).gt.icc)THEN WRITE(*,*)'error, icc is small!!' STOP ENDIF IR(j)=IR(j)+1 IF(IR(j).gt.icc)THEN WRITE(*,*)'error, icc is small!!' STOP ENDIF CM(i,IR(i))=sx CM(j,IR(j))=sx ICM(i,IR(i))=j ICM(j,IR(j))=i ENDIF ENDDO ENDDO DO i=1,m DIA(i)=omega/CM(i,1) SN(i)=IR(i) ENDDO kirmax=SN(idamax(m,SN,1)) omega1=1.-omega gm1=1./gamma icont=0 nacc=0 ! !********Start of outer method********** ! ! ! Compute (A-gamma D) x + q ! DO i=1,n HX(i)=Q(i) ENDDO DO j=2,iban-2,3 DO i=1,n HX(i)=HX(i)+A(i,j)*X(IA(i,j))+A(i,j+1)*X(IA(i,j+1))+A(i,j+2)* 1 X(IA(i,j+2)) ENDDO ENDDO DO j=2+((iban-1)/3)*3,iban DO i=1,n HX(i)=HX(i)+A(i,j)*X(IA(i,j)) ENDDO ENDDO DO i=1,n HXOLD(i)=HX(i)+A(i,1) * X(i) HX(i)=HX(i)+(A(i,1)-gamma*DD(i)) * X(i) ENDDO iter=1 ! ! outer iteration ! 999 continue DO i=1,n S1(i)=HX(i) * gm1 *AINV(i) ENDDO DO i=1,m S(i)=-B(i) ENDDO DO j=1,ibanc-2,3 DO i=1,m S(i)=S(i)-C(i,j)*S1(IC(i,j))-C(i,j+1)*S1(IC(i,j+1))- 1 C(i,j+2)*S1(IC(i,j+2)) ENDDO ENDDO DO j=1+(ibanc/3)*3, ibanc DO i=1,m S(i)=S(i)-C(i,j)*S1(IC(i,j)) ENDDO ENDDO DO i=1,m S(i)=S(i)*gamma ENDDO ! ! projected SOR ! DO 777 iter1=1,maxit1 DO j=1,m YZ0(j)=YZ(j) ENDDO DO i=1,m sum=-S(i) DO j=2,ir(i) sum=sum-CM(i,j)*YZ(ICM(i,j)) ENDDO YZ(i)=YZ0(i)*omega1+sum*DIA(i) if(i.le.mdq) YZ(i)=max(0.d0,YZ(i)) SN(i)=(-sum+CM(i,1)*YZ(i))*gm1 ENDDO ryz=0.0 DO i=1,mdq ryz=ryz+SN(i)*YZ(i) ENDDO DO i=1,m YZ0(i)=YZ(i)-YZ0(i) ENDDO IF((dnrm2(m,YZ0,1).lt.epsk/10*dnrm2(m,YZ,1).and.abs(ryz).lt.epsk) + .and. iter1.gt.2) THEN icont=icont+iter1 goto 800 ENDIF 777 continue info=2 icont=icont+maxit1 800 continue DO i=1,n TEMP(i)=0. ENDDO DO j=1,m DO i=1,kc(j) TEMP(IC(j,i))=TEMP(IC(j,i))+C(j,i)*YZ(j) ENDDO ENDDO DO i=1,n S1(i)=TEMP(i)-HX(i) ENDDO DO i=1,n X(i)=S1(i)*gm1*AINV(i) ENDDO DO i=1,n S1(i)=0.0 TTT(i)=X(i)-X0(i) ENDDO DO j=1,iban DO i=1,n S1(i)=S1(i)+A(i,j)*TTT(IA(i,j)) ENDDO ENDDO ! ! check if || A*d||eps*dnrm2(n,TTT(1),1)**2) then ! ! update of gamma ! gamma=0.0 xnum=0.0 IF(indrule.eq.1)THEN DO i=1,n gamma=gamma+S1(i)*AINV(i)*S1(i) xnum=xnum+TTT(i)*S1(i) ENDDO gamma=gamma/xnum ELSE DO i=1,n gamma=gamma+TTT(i)*AINV(i)*TTT(i) xnum=xnum+TTT(i)*S1(i) ENDDO gamma=xnum/gamma ENDIF c c compute theta_k c teta=0.0 DO i=1,n teta=teta+ X0(i)*S1(i) ENDDO DO i=1,n teta=teta+Q(i)*TTT(i) ENDDO teta=-teta/xnum IF(teta.lt.1 .and. teta.gt.0d0)THEN nacc=nacc+1 DO i=1,n X(i)=X0(i)+teta*TTT(i) ENDDO ELSE teta=1.0 ENDIF gm1=1./gamma ELSE teta=1 ENDIF DO i=1,n HXOLD(i)=HXOLD(i)+teta*S1(i) HX(i)=HXOLD(i)-gamma*X(i)*DD(i) SN(i)=HXOLD(i)-TEMP(i) ENDDO e3=abs(SN(idamax(n,SN,1))) epsk=min(epsk,max(e3**1.5,epssor)) e3=e3/rele3 DO i=1,n X0(i)=X0(i)-X(i) ENDDO ! ! outer stopping rule ! IF(e3.lt.eps.and.dnrm2(n,X0,1)/dnrm2(n,X(1),1).lt.eps)RETURN DO i=1,n X0(i)=X(i) ENDDO IF(iter .eq. maxit) THEN ! ! maximum number of outer interations ! info=1 RETURN ENDIF ! ! repeat the iteration ! iter=iter+1 GOTO 999 END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! PCG method with Arithmetic Mean Preconditioner ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE cg(CM,imm,ICM,IR,it,S,m,YZ, 1 maxit1,epsk,icc,SN,IK,D) PARAMETER(im=5001) IMPLICIT DOUBLE PRECISION(a-h,o-z) DIMENSION CM(imm,icc),ICM(imm,icc),IR(imm),S(imm),YZ(imm), 1 YZ0(imm),SN(imm),D(imm),IK(imm) DIMENSION P(im),AUX(im),T1(im),T2(im) ! DO i=1,m SN(i)=-S(i) ENDDO DO i=1,m DO j=1,IR(i) SN(i)=SN(i)-CM(i,j)*YZ(ICM(i,j)) ENDDO ENDDO ! ! preconditioner ! DO i=1,m T1(i)=SN(i) T2(i)=SN(i) ENDDO CALL lower(CM,ICM,imm,IK,T1,icc,m,D) CALL upper(CM,ICM,imm,IK,IR,T2,icc,m,D) DO i=1,m P(i)=0.5*(T1(i)+T2(i)) ENDDO ! xnum=ddot(m,P,1,SN,1) DO it=1,maxit1 DO i=1,m AUX(i)=0 DO j=1,IR(i) AUX(i)=AUX(i)+CM(i,j)*P(ICM(i,j)) ENDDO ENDDO xden=ddot(m,P,1,AUX,1) alfa=xnum/xden DO i=1,m YZ(i)=YZ(i)+alfa*P(i) SN(i)=SN(i)-alfa*AUX(i) T1(i)=SN(i) T2(i)=SN(i) ENDDO ryz=dabs(sn(idamax(m,SN,1))) IF(ryz.lt.epsk.and.it.gt.2)goto 800 CALL lower(CM,ICM,imm,IK,T1,icc,m,D) CALL upper(CM,ICM,imm,IK,IR,T2,icc,m,D) DO i=1,m AUX(i)=0.5*(T1(i)+T2(i)) ENDDO xden=xnum xnum=ddot(m,AUX,1,SN,1) beta=xnum/xden DO i=1,m P(i)=AUX(i)+beta*P(i) ENDDO ENDDO 800 continue RETURN END SUBROUTINE lower(CM,ICM,im,IK,X,icc,m,D) IMPLICIT DOUBLE PRECISION (a-h,o-z) INTEGER ICM(im,icc),IK(im) DIMENSION CM(im,icc),X(im),D(im) X(1)=X(1)*D(1) DO i=2,m DO j=2,IK(i) X(i)=X(i)-CM(i,j)*X(ICM(i,j)) ENDDO X(i)=X(i)*D(i) ENDDO RETURN END SUBROUTINE upper(CM,ICM,im,IK,IR,X,icc,m,D) IMPLICIT DOUBLE PRECISION (a-h,o-z) INTEGER ICM(im,icc),IK(im),IR(im) DIMENSION CM(im,icc),X(im),D(im) X(m)=X(m)*D(m) DO i=m-1,1,-1 DO j=IK(i)+1,IR(i) X(i)=X(i)-CM(i,j)*X(ICM(i,j)) ENDDO X(i)=X(i)*D(i) ENDDO RETURN END