PROGRAM NAG_TEST ! 24/05/98 ! This program computes by the NAG routine E04NKF 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_| ! ! In input, the user is only asked to introduce the name of the files ! containing the data. ! ! The program uses double precision floating point arithmetic. ! ! For some basic linear algebra operation, the program uses the ! BLAS1 routines. ! !--------------- DESCRIPTION OF THE INPUT VARIABLES ---------------------------------- ! ! name = name for the file containing the matrix A (character) ! namev = name for the file containing the constraints (character) ! ! 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 "nag.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 E04NKF and by the program. ! In particular, the routine E04NKf gives the ! following results: ! ! ifail = error flag ! X = vector of the computed solution; ! CLAMDA = 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 error 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 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit double precision(a-h,o-z) PARAMETER(in=5001,im=2001,ibmax=300,icmax=200, 1inm1=icmax*im+in,in1=in+1, 1inm=icmax*im,inq=ibmax*in,leniz=200000,lenz=13000000) ! ! 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 ! real secnds real*8 bl(IN+IM+1),bu(IN+IM+1),clamda(IN+IM+1),z(lenz),s1(IN) integer istate(IN+IM+1),iz(leniz) character*8 crname(IM+IN+1),names(5) COMMON/DATA/ C(IM,icmax),B(IM),Q(IN),A(IN,ibmax),XE(IN), + cc(IM*icmax+IN),x(IN+IM+1),indkc(IN+1),ihc(IM*icmax+IN), 1 ia(in,ibmax),ic(im,icmax),k1(in),k2(in),kc(im) EXTERNAL e04nkf,qphx data c/inm*0.0/,a/inq*0.0/,ic/inm*0/,ia/inq*0/ data b/im*0.0/,q/in*0.0/,cc/inm1*0.0/,ihc/inm1*0/,indkc/in1*0/ ! ! Executable statement ! call leggif(n,iban,m,meq,nva,ibanc,fmin,nrango,nnz, 1 sparsa,sparsc) m=m+1 do i=1,n bl(i)=-1d+25 bu(i)=1d+25 enddo bl(n+m)=-1d+25 bu(n+m)=1d+25 do i=n+1,n+m-1 bl(i)=b(i-n) enddo do i=n+1,n+m-1-meq bu(i)=1d+25 enddo do i=n+m-meq,n+m-1 bu(i)=bl(i) enddo do i=1,5 names(i)=' ' enddo do i=1,m+n crname(i)=' ' enddo do i=1,n istate(i)=0.0 x(i)=0.0 enddo open(unit=1,file='nag.dat',form='formatted') WRITE(1,*) WRITE(1,*) 'Results of E04NKF routine' WRITE(1,*) WRITE(1,*)'***INPUT PARAMETERS***' WRITE(1,*)'size of the problem =',n WRITE(1,*)'number of the linear constraints =', m-1 WRITE(1,*)'number of the equality constraints=',meq WRITE(1,*)'rank of A=',nrango 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,*) ifail=-1 call e04nmf('Print Level = 0') t1=secnds(0.0) call e04nkf(n,m,nnz,m,n,qphx,cc,ihc,indkc,bl,bu,'C', +names,n+m,crname,ns,x,istate,miniz,minz,ninf,sinf,obj, +clamda,iz,leniz,z,lenz,ifail) t2=secnds(real(t1)) write(1,*)'miniz=',miniz write(1,*)'minz=',minz write(1,*)'ifail=',ifail ! ! 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 WRITE(1,*)'elapsed time=',t2 ! ! Compute the active constraints in the solution ! DO i=1,m-meq-1 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-1 IF(ABS(XE(i)) .le. 1d-11) 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-n,clamda(i),i=n+1,n+m-1) WRITE(1,*) '***********************' STOP END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Matrix-vector product; matrix a is stored in compressed form ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine qphx(nstate,ncolh,xx,hx) implicit double precision(a-h,o-z) PARAMETER(IN=5001,IM=2001,ibmax=300,icmax=200, 1inm=icmax*im,inq=ibmax*in,leniz=200000,lenz=13000000) dimension xx(ncolh),hx(ncolh) COMMON/DATA/ C(IM,icmax),B(IM),Q(IN),A(IN,ibmax),XE(IN), +cc(IM*icmax+IN),x(IN+IM+1),indkc(IN+1),ihc(IM*icmax+IN), 1ia(in,ibmax),ic(im,icmax),k1(in),k2(in),kc(im) do i=1,ncolh hx(i)=0.0 enddo do i=1,ncolh do j=1,k1(i)+k2(i) hx(i)=hx(i)+a(i,j)*xx(ia(i,j)) enddo enddo return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Subroutine for reading the input data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine LEGGIF(N,iban,mm,meq,nva,ibanc,fmin,nrango,ind, 1 sparsa,sparsc) PARAMETER(IN=5001,IM=2001,ibmax=300,icmax=200) implicit double precision(a-h,o-z) COMMON/DATA/ C(IM,icmax),B(IM),Q(IN),A(IN,ibmax),ak(IN), +cc(IM*icmax+IN),x(IN+IM+1),indkc(IN+1),ihc(IM*icmax+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(5,'(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 1 i=1,n ia(i,1)=i k2(i)=0 1 k1(i)=1 do 5 i=1,inc/4 read(8) (l(j),m(j),coeff(j),j=1,4) do 10 j=1,4 if(l(j).eq.m(j))then a(l(j),1)=coeff(j) else k1(l(j))=k1(l(j))+1 a(l(j),k1(l(j)))=coeff(j) ia(l(j),k1(l(j)))=m(j) endif 10 continue 5 continue irest=mod(inc,4) if (irest. ne. 0) then read(8)(l(j),m(j),coeff(j),j=1,irest) do 15 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 15 continue endif read(8) (q(i),i=1,n) read(8) fmin read(8) sparsa close(8) do 3456 i=1,n do 3456 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 3456 continue do 77 i=1,n 77 ak(i)=k1(i)+k2(i) 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 11 i=1,mm 11 kc(i)=0 do 51 i=1,inc/4 read(8) (l(j),m(j),coeff(j),j=1,4) do 51 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) 51 continue irest=mod(inc,4) if (irest .ne. 0) then read(8)(l(j),m(j),coeff(j),j=1,irest) do 151 j=1,irest kc(l(j))=kc(l(j))+1 c(l(j),kc(l(j)))=coeff(j) ic(l(j),kc(l(j)))=m(j) 151 continue endif do i=1,mm b(i)=kc(i) enddo ibanc=b(idamax(mm,b,1)) ind=1 do j=1,n indkc(j)=ind do i=1,mm do jl=1,kc(i) if (Ic(i,jl).eq.j) then cc(ind)=c(i,jl) ihc(ind)=i ind=ind+1 goto 112 endif enddo 112 continue enddo cc(ind)=q(j) ihc(ind)=mm+1 ind=ind+1 enddo ind=ind-1 indkc(n+1)=ind+1 read(8) (ak(j),j=1,n) read(8) (b(j),j=1,mm) read(8) sparsc RETURN end