!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! PROGRAM AVPM (Adaptive 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 AVPM 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 AVPM are solved by a projected SOR ! method (Cryer, SIAM J. Control, 1971). ! In presence of equality constraints only, AVPM 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 inputavpm ! 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 AVPM ! (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) ! beta = parameter of the AVPM; beta must be a positive real ! less than 1; a suggested value for beta is 0.5; ! eta = parameter of the AVPM; eta must be a positive real ! greater than 0.5; a suggested value for eta is 0.9. ! ! 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 "avpm.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 AVPMSOR or AVPMPCG and by the program. ! In particular, the routine AVPMSOR or AVPMPCG gives the ! following results: ! ! info = if info=1, AVPM 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 AVPM iterations ! icont = total number of inner iterations ! nacc = additional projection steps performed in the search procedure of ! the projection parameter ! 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 AVPM 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 (AVPM) ! 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 READ(*,*) beta, eta ! ! 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 the matrix D ! 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(idamax(n,DD,1)) open(unit=1,file='avpm.dat',form='formatted') WRITE(1,*) WRITE(1,*) 'Adaptive 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 B =',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,*)'beta=', beta,' eta=', eta WRITE(1,*) IF(m-meq.ne.0)then t1=secnds(real(0.)) CALL AVPMSOR(n,m,meq,eps,maxit,maxit1,omega,info,iter, 1icont,iban,ibanc,kirmax,epssor,epstart,rele3,gamma, 2nacc,indrule,beta,eta) t2=secnds(real(t1)) ELSE t1=secnds(real(0.)) CALL AVPMPCG(n,m,meq,eps,maxit,maxit1,omega,info,iter, 1icont,iban,ibanc,kirmax,epssor,epstart,rele3,gamma, 2nacc,indrule,beta,eta) t2=secnds(real(t1)) ENDIF IF(info.eq.1) then WRITE(1,*) 'AVPM does not converge' STOP ENDIF WRITE(1,*)'info =',info WRITE(1,*)'Elapsed time in seconds =',t2 WRITE(1,*)'Number of AVPM 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 additional projection steps performed' WRITE(1,*)'in the search procedure of the correction', 1 ' parameter=',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 AVPM combined with the projected SOR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE AVPMSOR(n,m,meq,eps,maxit,maxit1,omega, 1info,iter,icont,iban,ibanc,kirmax,epssor,epstart, 2rele3,gamma,nacc,indrule,beta,eta) ! 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 ! a1=0.0d0 xnum=0.0d0 DO i=1,n a1=a1+X0(i)*S1(i)+Q(i)*TTT(i) xnum=xnum+TTT(i)*S1(i) ENDDO IF(a1.lt.0.and.-a1.lt.eta*xnum.and.iter.gt.1)then gamma=gamma/beta gm1=1./gamma DO i=1,n HX(i)=HXOLD(i)-gamma*DD(i)*X0(i) ENDDO nacc=nacc+1 goto 999 ELSE gamma=0.0 IF(indrule.eq.1)THEN DO i=1,n gamma=gamma+S1(i)*AINV(I)*S1(i) ENDDO gamma=gamma/xnum ELSE DO i=1,n gamma=gamma+TTT(i)*DD(i)*TTT(i) ENDDO gamma=eta*xnum/gamma ENDIF gm1=1.0d0/gamma ENDIF ENDIF DO i=1,n HXOLD(i)=HXOLD(i)+S1(i) HX(i)=HXOLD(i)-gamma*DD(i)*X(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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Subroutine of AVPM combined with the PCG method !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE AVPMPCG(n,m,meq,eps,maxit,maxit1,omega, 1info,iter,icont,iban,ibanc,kirmax,epssor,epstart, 2rele3,gamma,nacc,indrule,beta,eta) ! 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),IK(im),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=epssor ! ! compute the complementarity matrix; this matrix is ! stored in crompressed form in CM (ICM and KC) ! DO i=1,m IK(i)=1 ENDDO 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 IK(j)=IK(j)+1 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 ! ! PCG algorithm ! CALL cg(CM,IM,ICM,IR,iter1,S,m,YZ,maxit1,epsk,icc,SN, 1 IK,DIA) IF(iter1.eq.maxit1)THEN info=2 ENDIF icont=icont+iter1 ! 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||