-C ********************************************************* -C * * -C * UPTARR: VERIFIED * -C * * -C * Subroutine UPTARR: updates the arrays A, IA, JA, AT, * -C * IAT, JAT when a slack variable is added. * -C * * -C * called by: ISLACK subroutines called: NONE * -C ********************************************************* -C - SUBROUTINE UPTARR(iblk,COEFSK,ROW,C,ZEROD) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE - CHARACTER*8 NAMERW, NAMECL - DOUBLE PRECISION A,AT,COEFSK,C(SIZEVEC),ZEROD,Y,X,SLACK -C - INTEGER NROWS, NCOLS, IA, JA, IAT, JAT, ROW, - 1 iblk,iaclof,RWIAT -C - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof -C - IF(IBLK.EQ.1) THEN -C /* Push data up if iblk=1*/ -C /* Do IA,JA,A */ - DO 10 I1=IA(IBLA(2))+1,2,-1 - JA(I1)=JA(I1-1) - A(I1)=A(I1-1) - 10 CONTINUE - DO 15 I1=IBLA(2)+1,2,-1 - IA(I1)=IA(I1-1)+1 - 15 CONTINUE -C /* Now do IAT, JAT, AT */ -C /* JAT()++, Find row,JAT(.++)-JAT(end++),iat(.)++, iblat same */ -C - DO 20 I1=1,IAT(IBLAT(2))-1 - JAT(I1)=JAT(I1)+1 - 20 CONTINUE - RWIAT=IBLAT(1)-1+ROW - DO 25 I1=IAT(IBLAT(2))+1,IAT(RWIAT)+1,-1 - JAT(I1)=JAT(I1-1) - AT(I1)=AT(I1-1) - 25 CONTINUE - JAT(IAT(RWIAT))=1 - AT(IAT(RWIAT))=COEFSK - DO 30 I1=RWIAT+1,IBLAT(2) - IAT(I1)=IAT(I1)+1 - 30 CONTINUE -C - IBLA(2)=IBLA(2)+1 - NIA=NIA+1 - NJA=NJA+1 - NJAT=NJAT+1 - JA(1)=ROW - A(1)=COEFSK - IACLOF=IACLOF+1 - IBLSZ(1,1)=IBLSZ(1,1)+1 - NCOLS=NCOLS+1 - DO 32 I=NCOLS,2,-1 -C /* Slack gets moved up also*/ - SLACK(I)=SLACK(I-1) - NAMECL(I)=NAMECL(I-1) - C(I)=C(I-1) - 32 CONTINUE - SLACK(1)=-1.0D0*Y(ROW) - C(1)=ZEROD - NAMECL(1)='SLACK ' -C /* For W, add to end of data*/ - ELSE IF(IBLK.GT.ISCEN+1.AND.IBLK.LE.2*ISCEN+1) THEN -C /* IA,JA,A */ - NJA=NJA+1 - JA(NJA)=ROW - A(NJA)=COEFSK - NIA=NIA+1 - IA(NIA+1)=IA(NIA)+1 - IBLA(IBLK+1)=IBLA(IBLK+1)+1 -C /* IAT,JAT,AT */ - RWIAT=IBLAT(IBLK)+ROW-1 - DO 35 I=IAT(IBLAT(IBLK+1)),IAT(RWIAT+1)+1, -1 - JAT(I)=JAT(I-1) - AT(I)=AT(I-1) - 35 CONTINUE - JAT(IAT(RWIAT+1))=IBLSZ(IBLK,1)+1 - AT(IAT(RWIAT+1))=COEFSK - DO 40 I=RWIAT+1,IBLAT(IBLK+1) - IAT(I)=IAT(I)+1 - 40 CONTINUE - NJAT=NJAT+1 - NCOLS=NCOLS+1 - NAMECL(NCOLS)='SLACK ' - C(NCOLS)=ZEROD - IBLSZ(IBLK,1)=IBLSZ(IBLK,1)+1 -C -C If needed, I can also update the starting vector y to |c| <- (|c|^2+1)^0.5 -C and |A^T.b| <- (|A^T.b|+b_i)^0.5, and then multiply b through again by this -C factor and put the result in y. I don't think it will make a lot of diffence -C - ELSE - WRITE(7,*) '____________________________' - WRITE(7,*) 'NOW IN T, NOT ADDING SLACK UNTIL W. IBLK=',IBLK - WRITE(7,*) 'ROW=',ROW - WRITE(7,*) '______This mesg should never appear___________' - ENDIF - RETURN - END ##### sed > BQ_Decomp/impdrt.f <<'#####' 's/^-//' -c -C ********************************************************** -C * Subroutine IMPDRT, finds the direction of improvement * -C * DY. (pHASE 1) * -C * called by: Main Program subroutines called:MULTVS * -C * * -C * IMPUT PARAMETERS: dy, v, w, slack, nrows, ncols, b * -C * OUTPUT PARAMETERS: dy * -C ********************************************************** -C -C - SUBROUTINE IMPDRT(DY,V,W,D,B, RCOND) - INCLUDE (sizes.h) -C - DOUBLE PRECISION V(SIZEVEC), W(SIZEVEC), D(SIZEVEC), B(SIZEVEC) - DOUBLE PRECISION DY(SIZEVEC), PROD1, PROD2, AUX, AT, A, RCOND(4) -C - INTEGER IA,JA,IAT,JAT,NROWS,NCOLS -C - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof -C -C - DO 700 I=1,NROWS - V(I)= 0.0D0 - 700 CONTINUE -C /* Perform same func as FINDV*/ -C /* Code adapted from ATYCSL*/ -C /* Form v=AD2.e */ - DO 705 I=IBLAT(1),IBLAT(2)-1 - I1=I-IBLAT(1)+1 - DO 704 J=IAT(I),IAT(I+1)-1 - V(I1)=V(I1)+(D(JAT(J))*D(JAT(J)))*AT(J) - 704 CONTINUE - 705 CONTINUE - J2=IBLSZ(1,2) - DO 710 I=2,ISCEN+1 - DO 709 J=IBLA(I),IBLA(I+1)-1 - J1=(J-IBLA(I)+1)+IACLOF - DO 708 K=IA(J),IA(J+1)-1 - V(J2+JA(K))=V(J2+JA(K))+A(K)*D(J1)*D(J1) - 708 CONTINUE - 709 CONTINUE - J2=J2+IBLSZ(I,2) - 710 CONTINUE -C /* Now do Wty. I3=COL_OFS */ - I1=IBLSZ(1,2) - I2=IBLSZ(1,1) - DO 720 I=ISCEN+2,2*ISCEN+1 - DO 719 J=IBLA(I),IBLA(I+1)-1 - I3=I2+(J-IBLA(I)+1) - DO 718 K=IA(J),IA(J+1)-1 - V(I1+JA(K))=V(I1+JA(K))+A(K)*D(I3)*D(I3) - 718 CONTINUE - 719 CONTINUE - I1=I1+IBLSZ(I,2) - I2=I2+IBLSZ(I,1) - 720 CONTINUE -C...................................... -C /* Solve Bw=v */ - CALL FINDDY(-1,V,W, RCOND) -c CALL FINDDY(1,V,W, RCOND) -C /* Find vt.w */ - PROD1= 0.0D0 - DO 730 I=1,NROWS - PROD1= PROD1+V(I)*W(I) - 730 CONTINUE -C /* Solve vt.u */ - PROD2= 0.0D0 - DO 740 I=1,NROWS - PROD2= PROD2+V(I)*DY(I) - 740 CONTINUE -C /* Find et.D2.e */ - AUX= 0.0D0 - DO 750 I=1, NCOLS - AUX= AUX+D(I)*D(I) - 750 CONTINUE -C /* Set dy2, b(nrows+1)=-M */ - DY(NROWS+1)=(B(NROWS+1)+PROD2)/(AUX-PROD1) - -C /* Find Dy1 */ - DO 760 I= 1, NROWS - DY(I)= DY(I) + W(I)*DY(NROWS+1) - 760 CONTINUE -C - RETURN - END ##### sed > BQ_Decomp/main.f <<'#####' 's/^-//' -C *************************************************************** -C * AIX Version 3/3/92 * -C * * -C * Subroutines called: * -C * READAT, STARTY, CSLACK, CHECKY, SETPH1, INIMTR, GETSTX * -C * OUTPRD, SPRSPK, IJBEGN, INIJ, IJEND, ORDRA1, FREEVA * -C * INAIJ5, INBI, SOLVE5, MULATY, STEPSZ, DUALOB, ISLACK * -C * UPDMTR, DROPRW PRMLOB, PREPRO, PRINT, SORTCL, ARWISE * -C * DUALIZ, switch * -C * * -C *************************************************************** -C - INCLUDE (sizes.h) - CHARACTER*3 BLANK1,RTYPE,trtype - CHARACTER*4 PROBLM - CHARACTER*8 BLANK2, NAMERW, NAMECL, FEASBL,tnamec,tnamer - CHARACTER*70 argv -C - DOUBLE PRECISION A,B,C,Y,X,GAMMA,SLACK,STEP,DV(SIZEVEC),ZEROD, - 1 BIGMIU,LST2, SLKMIN, DY(SIZEVEC), ATY(SIZEVEC), DUMMY, - 2 ZOLD, YOLD, EPSILO, VALFIX, VALLBD, AAT, UBOUND, - 3 DIAG,LSTPRD,ZCONST, ZLOWER,ZUPPER, STOPC, RCOND(4), - 4 D(SIZEVEC),EPSUPD,UBSTEP,ARTFCL,BIGM, W(SIZEVEC), - 5 V(SIZEVEC), AT, ta, tb, tc, taux, EPSPH1 , AUX -C - REAL PREPTM, READTM, STUPTM, PHS1TM, PHS2TM -C - INTEGER IA,JA,NROWS,NCOLS,IAAT,JAAT,PTRPRD, ITERA1, ITERA2 - INTEGER ZEROI, LSTDST,MSGLVA,IERRA,MAXSA,NEQNS, PHASE, - 1 VARFIX, VARLBD, VRFREE, ORGIDX, NFREE, NFIX, NSTFRE, - 1 NLBD, NSTVAR, NRWDRP, VRCSGN, NCHSGN, ITERTN, - 1 IAT, JAT, ibla, iblat, iblaat, idiag, iscen, narows, - 1 nwrows, nwcols, itotc, itotr, tsort, nsort,toridx - INTEGER ERROR, tnrows, tncols, tia, tja, PTR2, DST2, JSTART -C - INTEGER I, itim(10) , ICPAD , IPRNTE, IPRNTS, - 1 IUNIT , MAPPTR, MAXINT, LEN, - 1 MXUSED, NVARS , STAGE, ISAVE, IDSAVE - REAL MCHEPS, RATIOL, RATIOS, SVALUE, CLOCK, RSAVE - DOUBLE PRECISION S(1), DSAVE, STAT(20) -C - COMMON /debug/ idbg - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /INPT2/ B(SIZEVEC), C(SIZEVEC), ZCONST - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - COMMON /OUTPRO/ LSTPRD(SIZEVEC),PTRPRD(SIZEMAT),LSTDST(SIZEMAT) - COMMON /OUTPR2/ IBLPTR(MAXBLK),IBLLST(MAXBLK),IBLDST(MAXBLK), - 1 LST2(SIZEPRD),DST2(SIZEPRD), PTR2(SIZEVEC) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /SPAUSR/ MSGLVA,IERRA,MAXSA,NEQNS - COMMON /TMPAAT/ NZTOTL,NIAAT - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, - 1 STUPTM, PHS1TM, PHS2TM, MXITNS - COMMON /CONSTR/ NAMERW(SIZEVEC),RTYPE(SIZEVEC),NAMECL(SIZEVEC) - COMMON /FIXED/ VALFIX(SIZEBLK),VALLBD(SIZEBLK),VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK),VRFREE(SIZEVEC),VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, - 1 MCHEPS, CLOCK - COMMON /SPACON/ STAGE , MXUSED, ICPAD(48) - COMMON /SPAMAP/ MAPPTR(50) - COMMON /SPADTA/ SVALUE(50) -C - COMMON /SPSAVE/ ISAVE(MAXBLK,103), RSAVE(MAXBLK,50), - 1 IDSAVE(MAXBLK), DSAVE(SIZEDSAVE) -C -C A(I) : contains the nonzero elements of the matrix A -C B : dense representation of the resource vector b. -C C : dense representation of the objective cost vector c. -C Y(I) : current solution vector -C X(I) : tentative solution of the primal problem -C D(I) : aproximates the diagonal els. of Dk=diag(1/SLACK(I)) -C DY : vecto DY, direction of improvement at each iteration -C IA(I) : contains pointers to the first entry of each column -C The component IA(n+1) points to the position after -C the last entry in A. -C JA(I) : contains th row index of the nonzero element in A(I) -C NROWS : the number of rows of the coefficient matrix A. -C NCOLS : the number of columns of the coefficient matrix A. -C SLACK(I): current slack at the ith constraint -C PREPTM: preprocessing time -C READTM: input data reading time -C STUPTM: set up time before starting the iterations -C PHS1TM: time spent during the iterations of phase 1 -C PHS2TM: time spent during the iterations of phase 2 -C ITERA1: number of iterations in phase1 -C ITERA2: number of iterations in phase 2 -C ZOLD : most recent lower bnd on the optimal objective value -C EPSILO: relative improvement of the objctive value -C EPSUPD: minimal relative change in D(I) to update it. -C UBOUND: upper bound on the objective function value. -C UBSTEP: upper bound on the step size. -C GAMMA : PARAMETER -C VARFIX: variables whose values are already known. -C VRFREE: stores the indices (original) of the free variables -C VARLBD: stores the indices (original) of LBounded variables. -C VRCSGN: stores indices (original) variab. which sign changed -C VALFIX: stores the values of the fixed variables. -C VALLBD: stores the lower bounds of the LB variables. -C NFIX : number of fixed variables. -C NLBD : number of lower bounded variables. -C NFREE : number of unrestricted variables. -C NSTFRE: number of initial free variables. -C NCHSGN: number of variables whose sign were changed. -C NRWDRP: number of rows dropped (due to redundancy, etc.). -C NSTVAR: number of starting variables. -C ORGIDX: stores the original index of each column. -C PROBLM: specifies the problem being solved by the algorithm. -C FEASBL: the probl. is FEASIBLE, UNBOUNDED OR INFEASIBLE. -C ARTFCL: value of the artificial variable. -C idbg : debug level 0 is lowest, 4 is highest -C istrtc: starting column for a block -C iendc : ending column for a block -C istrtr: starting row for a block -C iendr : ending row for a block -C iscen : number of scenarios -C itotc : total number of columns read thus far from tia(.) -C itotr : total number of rows read thus far -C itotel: total number of array elements read thus far. -C narows: number of rows in A matrix -C nwrows: number of rows in W matrix -C nacols: number of columns in A and T matrix -C nwcols: number of columns in W matrix -C tia : temporary ia matrix: used for initial input -C tja : temporary ja matrix: used for initial input -C ta : temporary a matrix: used for initial input -C nia : number of elements in ia thus far. -C irwofs: row offset: used to interface sparse representations -C iaclof: column offset: for T matrices -C -C The COMMON block /SPAUSR/ allows this program to interact -C with the computer package SPARSPAK. Some subroutines of this -C package are used for solving the large scale linear system: -C (ADkDkAt)dy = b, at each iteration. The variables in this -C block are specified in the SPARSPAK manual. -C -C -------------------------- -C Initialize SPARSPAK-A -C -------------------------- - ITIM(1)=MCLOCK() - CALL SPRSPK - MAXSA=SIZEVEC - MSGLVA=0 - idbg=0 -C --------------------- -C INITIALIZE PARAMETERS -C --------------------- - i=dfhtime(1) - CALL GETARG(1,ARGV) - open(file=argv,unit=4,status='unknown') - CALL INTLZE(EPSILO,EPSUPD,ERROR,FEASBL,GAMMA,NFIX,NFREE, - 1 NLBD,NRWDRP,NSTVAR,PHASE,PROBLM,UBOUND, - 1 UBSTEP,YOLD,ZCONST,ZOLD,NCHSGN,STOPC,BIGMIU, EPSPH1) -C - CALL READAT(JSTART,EPSPH1,EPSILO,BIGMIU) - ITIM(2)=MCLOCK() - BIGMIU=1.D05 - JSTART=0 - call getarg(7,argv) - read(argv,*) BIGMIU - call getarg(8,argv) - READ(argv,*) JSTART -C - write(0,*) ' narows, nacols, nwrows, nwcols', - 1 narows, nacols, nwrows, nwcols - IAAT(1)=1 - NSTFRE= NFREE - CALL SORTCL - CALL STARTY(JSTART) -C . utput densities -c Do 1499 i1=1,nacols -c i3=0 -c do 1498 i2=tia(i1),Tia(i1+1)-1 -c if(tja(i2).ge.NAROWS+1.and.dabs(ta(i2)).ge.1.D-06) i3=i3+1 -c1498 continue -c write(6,*) ' COL: ',I1,' : ',TNAMEC(i1),' : ',i3 -c1499 continue -C -c write(6,*) ' ISCEN, NWCOLS, etc.',ISCEN,NACOLS,NAROWS,NWCOLS -C............................................. BEGIN FORVECION LOOP .......... -C - DO 15 iblk=1,2*ISCEN+1 -C - CALL SWITCH(iblk) -C - if(iblk.eq.1) then - irwofs=zeroi - iclofs=zeroi - else if(iblk.le.iscen+1) then - irwofs=0 - do 1 i=1,iblk-1 - irwofs=irwofs+iblsz(i,2) - 1 continue - iclofs=iaclof - else - irwofs=iblsz(1,2) - iclofs=iblsz(1,1) - do 2 i=iscen+2,iblk-1 - irwofs=irwofs+iblsz(i,2) - iclofs=iclofs+iblsz(i,1) - 2 continue - endif -c - IF (FEASBL .NE. 'FEASIBLE') stop - CALL FREEVA(iblk,iclofs,irwofs,PROBLM,NFREE,VRFREE,ORGIDX,C) - CALL ISLACK(iblk,irwofs,PROBLM,VRFREE,ORGIDX,NFREE, - 1 C,ZEROD,NROWS) -C -C Since finding a P1 solution doesn't require modifying the A or AT data -C structures, I'm moving CSLACK and SETPH1 outside the loop that builds -C A,AT,AAT, etc. THen, the subs will be modified to be consistent with -C the new data structure. -C - 5 CONTINUE - IF(IBLK.GE.2.AND.IBLK.LE.ISCEN+1) THEN - IBLAAT(IBLK+1)=NIAAT - IDIAG(IBLK+1)=IDIAG(IBLK) - IBLPTR(IBLK+1)=IBLPTR(IBLK) - IBLLST(IBLK+1)=IBLLST(IBLK) - IBLDST(IBLK+1)=IBLDST(IBLK) - GOTO 7 - ENDIF - CALL STJAAT(iblk,ZEROI) - CALL OUTPRD(iblk,ERROR) - IF (ERROR .GT. ZEROI) then - write(0,*) 'after OUTPRD, error = ',error - stop - endif -C - CALL BLDMTR(iblk,iclofs,IRWOFS,D,ZEROD,EPSILO) -C - 7 CONTINUE -C -C - 15 CONTINUE - ITIM(3)=MCLOCK() -C -C -C ........................................ END FORMATION LOOP ............. - CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) - IF (SLKMIN .GT. ZEROD) GOTO 42 - PHASE=1 - CALL SETPH1(B,Y,SLKMIN,NROWS,BIGM,BIGMIU) - CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) -C /* If in P1, need to recalc D*/ -C /* Since BLDMTR used old slacks, update*/ - CALL UPDMTR(D,SLACK,-1.0D0) - WRITE(6,*) ' PHASE 1 NECESSARY, MIN SLACK=',SLKMIN - 42 CONTINUE -C ........................................ BEGIN SOLUTION LOOP ............ -C - i=dfhtime(1) - ITIM(9)=0 -C - DO 100 ITERTN=1,50 -c write(7,*) ' ............. SLACK ............. ' -c write(7,'(5E13.6)') (slack(i),i=1,NROWS+1) - -C -c write(6,*) ' ITER: ',ITERTN, ' yart: ',Y(NROWS+1) - do 880 i1=1,15 - stat(i1)=0.0D0 - 880 continue - i2=iblsz(1,1) - stat(3)=1.D+30 - stat(4)=1.D-30 - do 890 i1=1,iblsz(1,1) - stat(1)=stat(1)+SLACK(i1) - stat(2)=stat(2)+SLACK(i1)*SLACK(i1) - if(stat(3).ge.slack(i1)) stat(3)=slack(i1) - if(stat(4).le.slack(i1)) stat(4)=slack(i1) - 890 continue - stat(7)=1.D+30 - stat(8)=1.D-30 - do 891 i1=iblsz(1,1)+1,ncols - stat(5)=stat(5)+SLACK(i1) - stat(6)=stat(6)+SLACK(i1)*SLACK(i1) - if(stat(7).ge.slack(i1)) stat(7)=slack(i1) - if(stat(8).le.slack(i1)) stat(8)=slack(i1) - 891 continue - stat(9)=stat(1)/iblsz(1,1) - stat(10)=stat(5)/(ncols-iblsz(1,1)) - stat(11)=(stat(1)+stat(5))/(ncols) - stat(12)=stat(2)+stat(6) - stat(13)=stat(3) - if(stat(7).le.stat(13)) stat(13)=stat(7) - stat(14)=stat(4) - if(stat(8).ge.stat(14)) stat(14)=stat(8) - write(6,892) itertn,stat(9),(stat(2)-stat(9)*stat(9))/(i2-1), - 1 stat(3),stat(4) - write(6,893) itertn,stat(10),(stat(6)-stat(10)*stat(10))/ - 1 (ncols-i2-1),stat(7),stat(8) - write(6,894) itertn,stat(11),(stat(12)-stat(11)*stat(11))/ - 1 (ncols-1), stat(13),stat(14) - 892 format('A :It:',I3,' AV: ',E12.5,' VA: ',E12.5,' Lo: ',E12.5, - 1 ' Hi: ',E12.5) - 893 format('W :It:',I3,' AV: ',E12.5,' VA: ',E12.5,' Lo: ',E12.5, - 1 ' Hi: ',E12.5) - 894 format('All:It:',I3,' AV: ',E12.5,' VA: ',E12.5,' Lo: ',E12.5, - 1 ' Hi: ',E12.5) -C - ITIM(7)=MCLOCK() - CALL FINDDY(ITERTN,B,DY,RCOND) - ITIM(8)=MCLOCK() -C - ITIM(9)=ITIM(9)+ITIM(8)-ITIM(7) -C /* Following: impdrt,findv*/ - IF(PHASE.EQ.1) THEN - CALL IMPDRT(DY,V,W,D,B, RCOND) - ENDIF -C /* Now form Atdy */ -C /* must transfer dt to tmp*/ -c write(7,*) ' ITER: ',ITERTN -c write(7,*) ' DY = ' -c write(7,'(5(E13.6))') (dy(i),i=1,NROWS+1) - CALL ATYCSL(2,DY,DV,PHASE,SLKMIN) -c write(7,*) ' DV = ' -c write(7,'(5(E13.6))') (dv(i),i=1,NCOLS+1) - -C - CALL STEPSZ(STEP,UBSTEP,GAMMA,SLACK,DV,NCOLS,STOPC) -c write(7,*) ' STEP: ', STEP -C - if(itertn.eq.1) write(0,*) ' .............', - 1 ' Iteration Log ..............' -c if(MOD(itertn,5).eq.0) then - if(STOPC.ge.99) write(0,1192) itertn, phase, RCOND(1), RCOND(2) - if(STOPC.lt.99) write(0,1191) itertn, stopc, - 1 phase, RCOND(1), RCOND(2) -c endif -C - IF ((STEP/GAMMA) .EQ. UBSTEP) THEN - FEASBL= 'STUCK-1 ' - WRITE(0,*) 'Stuck1:(Step/Gamma),UBStep',STEP/GAMMA,UBSTEP - GOTO 120 - END IF -C - DO 80 I= 1, NROWS - Y(I)= Y(I) + STEP*DY(I) - 80 CONTINUE - IF (PHASE .EQ. 2) GOTO 90 - Y(NROWS+1)= Y(NROWS+1) + STEP*DY(NROWS+1) - write(0,*) ' yart: ',y(nrows+1),yold,yold-y(nrows+1) - IF (Y(NROWS+1) .GT. EPSILO) GOTO 85 - IF (Y(NROWS+1) .GT. 0.0D0) GOTO 90 -C - ITERA1= ITERTN - ITIM(4)=MCLOCK() - PHASE=2 -C - GOTO 90 -C - 85 IF ((YOLD - Y(NROWS+1)) .GE. EPSPH1) THEN - YOLD=Y(NROWS+1) - GOTO 95 - ENDIF - FEASBL= 'STUCK-2 ' - WRITE(0,*) 'Stuck2:(delta(y))=',YOLD - Y(NROWS+1) -C /* Calc DUALOB */ - ZLOWER= 0.0D0 - DO 86 I=1,NROWS - ZLOWER= ZLOWER + B(I) * Y(I) - 86 CONTINUE - GOTO 120 -C /* Calc DUALOB */ - 90 ZLOWER= 0.0D0 - DO 91 I=1,NROWS - ZLOWER= ZLOWER + B(I) * Y(I) - 91 CONTINUE - IF (ZLOWER .GT. UBOUND) THEN - FEASBL= 'UNBOUNDD' - WRITE(0,*) ' Unbounded' - GOTO 120 - ENDIF -c -c A note regarding stopping criterion: Here, we use the criterion -c that the obj function does not change in the 8th sig. digit. -c To do comparisons w/ OSL, need to decrease RMUlimit to 0. -c Note: Max. comp. slackness doesn't work (alg. not numerically -c stable enough.) (i.e. max_{i} (x(i)*slack(i)) < stpeps -c - AUX=DABS(ZOLD) - IF (AUX .LT. 1.0D0) AUX=1.0D0 - STOPC=DABS(ZLOWER-ZOLD)/AUX - IF (STOPC .LE. EPSSTP) GOTO 120 - ZOLD=ZLOWER -C /* Find new slacks */ - 95 CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) -C - CALL UPDMTR(D,SLACK,EPSUPD) - - 100 CONTINUE -C - 120 ITIM(5)=MCLOCK() - i=dfhtime(2) - IF (PHASE .EQ. 2) THEN - ITERA2= ITERTN-ITERA1 - ELSE - ARTFCL= Y(NROWS+1) - ITERA1= ITERTN - ENDIF -C /* Find Primal Objective */ - ZEROD =0.0D0 - ZUPPER=ZEROD - DO 122 I=1,NCOLS - X(I)= DV(I)/(SLACK(I)*SLACK(I)) - ZUPPER= ZUPPER+C(I)*X(I) - 122 CONTINUE -C /* GetZ value */ - IF (PROBLM .EQ. 'DUAL') GOTO 125 - DUMMY= ZUPPER - ZUPPER= -ZLOWER - ZLOWER= -DUMMY -C - 125 ZUPPER= ZUPPER + ZCONST - ZLOWER= ZLOWER + ZCONST -C - CALL GETSTX(X,NCOLS,ZEROI,NROWS,Y,PROBLM) -C 160 CALL PRINT(ZUPPER,ZLOWER,FEASBL,PROBLM,NFIX,NRWDRP, -C 1 NSTFRE,PHASE,ARTFCL,BIGM) -C - ITIM(6)=MCLOCK() - write(0,167) NZTOTL, ZLOWER - write(0,168) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), - 1 itim(5)-itim(4),itim(6)-itim(5),itim(9) - write(0,169) itim(6)-itim(1) - write(6,167) NZTOTL, ZLOWER - write(6,168) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), - 1 itim(5)-itim(4),itim(6)-itim(5),itim(9) - write(6,169) itim(6)-itim(1) - write(6,170) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), - 1 itim(5)-itim(4),itim(6)-itim(5),itim(9),itim(6)-itim(1),NIAAT, - 2 ZLOWER - - 167 format(/' Nzs in AAT :',i15,' Solution :',E15.9) - 168 format(' Reading Time :',i15,' Formation Time :',i15/ - 1 ' Phase I Time :',i15,' Phase II Time :',i15/ - 2 ' Print Time :',i15,' TOTAL finddy time :',i15) - 169 format(' TOTAL TIME: : ',i15) - 170 format('SUM: ',6(I12,','),',,',2I12,E15.9) - 1191 format('Iteration: ',i5,' Stop Criterion: ',F15.12, - 1 ' Phase:',I2,' kap(G1): ',E18.10,' kap(G2): ',E18.10) - 1192 format('Iteration: ',i5,' Phase:',i2, ' kap(G1): ',E18.10, - 1 ' kap(G2): ',E18.10) - 180 STOP - END ##### sed > BQ_Decomp/mainbak.f <<'#####' 's/^-//' -C -C *************************************************************** -C * * -C * AIX Version 3/3/92 -C * * -C * This program * -C * (1) Reads in the coefficients of a LP program in the * -C * standard form min cx, s.t. Ax=b, x>=0. * -C * (2) Solves the LP program using the DUAL AFFINE * -C * scaling version of the Karmarkar's algorithmn * -C * (3) Prints out the final solution * -C * * -C * It uses, at each iteration, some subroutines of the * -C * SPARSPAK-A software package to find the direction of * -C * improvement. * -C * * -C * Subroutines called: * -C * READAT, STARTY(JSTART), CSLACK, CHECKY, SETPH1, INIMTR, GETSTX * -C * OUTPRD, SPRSPK, IJBEGN, INIJ, IJEND, ORDRA1, FREEVA * -C * INAIJ5, INBI, SOLVE5, MULATY, STEPSZ, DUALOB, ISLACK * -C * UPDMTR, DROPRW PRMLOB, PREPRO, PRINT, SORTCL, ARWISE * -C * DUALIZ, switch * -C * * -C *************************************************************** -C - INCLUDE (sizes.h) - CHARACTER*3 BLANK1,RTYPE,trtype - CHARACTER*4 PROBLM - CHARACTER*8 BLANK2, NAMERW, NAMECL, FEASBL,tnamec,tnamer - CHARACTER*70 argv -C - DOUBLE PRECISION A,B,C,Y,X,GAMMA,SLACK,STEP,DV(SIZEVEC),ZEROD, - 1 BIGMIU,LST2, SLKMIN, DY(SIZEVEC), ATY(SIZEVEC), DUMMY, - 2 ZOLD, YOLD, EPSILO, VALFIX, VALLBD, AAT, UBOUND, - 3 DIAG,LSTPRD,ZCONST, ZLOWER,ZUPPER, STOPC, - 4 D(SIZEVEC),EPSUPD,UBSTEP,ARTFCL,BIGM, W(SIZEVEC), - 5 V(SIZEVEC), AT, ta, tb, tc, taux, EPSPH1 , aux2 -C - REAL PREPTM, READTM, STUPTM, PHS1TM, PHS2TM -C - INTEGER IA,JA,NROWS,NCOLS,IAAT,JAAT,PTRPRD, ITERA1, ITERA2 - INTEGER ZEROI, LSTDST,MSGLVA,IERRA,MAXSA,NEQNS, PHASE, - 1 VARFIX, VARLBD, VRFREE, ORGIDX, NFREE, NFIX, NSTFRE, - 1 NLBD, NSTVAR, NRWDRP, VRCSGN, NCHSGN, ITERTN, - 1 IAT, JAT, ibla, iblat, iblaat, idiag, iscen, narows, - 1 nwrows, nwcols, itotc, itotr, tsort, nsort,toridx - INTEGER ERROR, tnrows, tncols, tia, tja, PTR2, DST2, JSTART -C - INTEGER I, itim(10) , ICPAD , IPRNTE, IPRNTS, - 1 IUNIT , MAPPTR, MAXINT, LEN, - 1 MXUSED, NVARS , STAGE, ISAVE, IDSAVE - REAL MCHEPS, RATIOL, RATIOS, SVALUE, CLOCK, RSAVE - DOUBLE PRECISION S(1), DSAVE -C - COMMON /debug/ idbg - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /INPT2/ B(SIZEVEC), C(SIZEVEC), ZCONST - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - COMMON /OUTPRO/ LSTPRD(SIZEVEC),PTRPRD(SIZEMAT),LSTDST(SIZEMAT) - COMMON /OUTPR2/ IBLPTR(MAXBLK),IBLLST(MAXBLK),IBLDST(MAXBLK), - 1 LST2(SIZEPRD),DST2(SIZEPRD), PTR2(SIZEVEC) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /SPAUSR/ MSGLVA,IERRA,MAXSA,NEQNS - COMMON /TMPAAT/ NZTOTL,NIAAT - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, - 1 STUPTM, PHS1TM, PHS2TM, MXITNS - COMMON /CONSTR/ NAMERW(SIZEVEC),RTYPE(SIZEVEC),NAMECL(SIZEVEC) - COMMON /FIXED/ VALFIX(SIZEBLK),VALLBD(SIZEBLK),VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK),VRFREE(SIZEVEC),VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, - 1 MCHEPS, CLOCK - COMMON /SPACON/ STAGE , MXUSED, ICPAD(48) - COMMON /SPAMAP/ MAPPTR(50) - COMMON /SPADTA/ SVALUE(50) -C - COMMON /SPSAVE/ ISAVE(MAXBLK,103), RSAVE(MAXBLK,50), - 1 IDSAVE(MAXBLK), DSAVE(SIZEDSAVE) -C -C A(I) : contains the nonzero elements of the matrix A -C B : dense representation of the resource vector b. -C C : dense representation of the objective cost vector c. -C Y(I) : current solution vector -C X(I) : tentative solution of the primal problem -C D(I) : aproximates the diagonal els. of Dk=diag(1/SLACK(I)) -C DY : vecto DY, direction of improvement at each iteration -C IA(I) : contains pointers to the first entry of each column -C The component IA(n+1) points to the position after -C the last entry in A. -C JA(I) : contains th row index of the nonzero element in A(I) -C NROWS : the number of rows of the coefficient matrix A. -C NCOLS : the number of columns of the coefficient matrix A. -C SLACK(I): current slack at the ith constraint -C PREPTM: preprocessing time -C READTM: input data reading time -C STUPTM: set up time before starting the iterations -C PHS1TM: time spent during the iterations of phase 1 -C PHS2TM: time spent during the iterations of phase 2 -C ITERA1: number of iterations in phase1 -C ITERA2: number of iterations in phase 2 -C ZOLD : most recent lower bnd on the optimal objective value -C EPSILO: relative improvement of the objctive value -C EPSUPD: minimal relative change in D(I) to update it. -C UBOUND: upper bound on the objective function value. -C UBSTEP: upper bound on the step size. -C GAMMA : PARAMETER -C VARFIX: variables whose values are already known. -C VRFREE: stores the indices (original) of the free variables -C VARLBD: stores the indices (original) of LBounded variables. -C VRCSGN: stores indices (original) variab. which sign changed -C VALFIX: stores the values of the fixed variables. -C VALLBD: stores the lower bounds of the LB variables. -C NFIX : number of fixed variables. -C NLBD : number of lower bounded variables. -C NFREE : number of unrestricted variables. -C NSTFRE: number of initial free variables. -C NCHSGN: number of variables whose sign were changed. -C NRWDRP: number of rows dropped (due to redundancy, etc.). -C NSTVAR: number of starting variables. -C ORGIDX: stores the original index of each column. -C PROBLM: specifies the problem being solved by the algorithm. -C FEASBL: the probl. is FEASIBLE, UNBOUNDED OR INFEASIBLE. -C ARTFCL: value of the artificial variable. -C idbg : debug level 0 is lowest, 4 is highest -C istrtc: starting column for a block -C iendc : ending column for a block -C istrtr: starting row for a block -C iendr : ending row for a block -C iscen : number of scenarios -C itotc : total number of columns read thus far from tia(.) -C itotr : total number of rows read thus far -C itotel: total number of array elements read thus far. -C narows: number of rows in A matrix -C nwrows: number of rows in W matrix -C nacols: number of columns in A and T matrix -C nwcols: number of columns in W matrix -C tia : temporary ia matrix: used for initial input -C tja : temporary ja matrix: used for initial input -C ta : temporary a matrix: used for initial input -C nia : number of elements in ia thus far. -C irwofs: row offset: used to interface sparse representations -C iaclof: column offset: for T matrices -C -C The COMMON block /SPAUSR/ allows this program to interact -C with the computer package SPARSPAK. Some subroutines of this -C package are used for solving the large scale linear system: -C (ADkDkAt)dy = b, at each iteration. The variables in this -C block are specified in the SPARSPAK manual. -C -C -------------------------- -C Initialize SPARSPAK-A -C -------------------------- - CALL SPRSPK - MAXSA=SIZEVEC - MSGLVA=0 - idbg=0 -C --------------------- -C INITIALIZE PARAMETERS -C --------------------- - CALL GETARG(1,ARGV) - write(6,*) ' argv(1)=',argv - open(file=argv,unit=4,status='unknown') - CALL INTLZE(EPSILO,EPSUPD,ERROR,FEASBL,GAMMA,NFIX,NFREE, - 1 NLBD,NRWDRP,NSTVAR,PHASE,PROBLM,UBOUND, - 1 UBSTEP,YOLD,ZCONST,ZOLD,NCHSGN,STOPC,BIGMIU, EPSPH1) -C - CALL READAT(JSTART,EPSPH1,EPSILO,BIGMIU) - BIGMIU=1.D05 - JSTART=0 -c call getarg(7,argv) -c read(argv,*) BIGMIU -c call getarg(8,argv) -c READ(argv,*) JSTART -C - IAAT(1)=1 -C - NSTFRE= NFREE - CALL SORTCL - CALL STARTY(JSTART) -C -c PREPTM= DTIME(0) - write(6,*) ' ISCEN, NWCOLS, etc.',ISCEN,NACOLS,NAROWS,NWCOLS - write(6,*) 'tCncols=',tncols -c do 191 i=1,tncols -c write(8,*) ' ia(',i,')=',tia(i) -c191 continue -c do 192 i=1,tia(tncols) -c write(8,*) 'ja(',i,')=',tja(i) -c192 continue -C............................................. BEGIN FORVECION LOOP .......... -C - DO 15 iblk=1,2*ISCEN+1 -C - CALL SWITCH(iblk) -C - if(iblk.eq.1) then - irwofs=zeroi - iclofs=zeroi - else if(iblk.le.iscen+1) then - irwofs=0 - do 1 i=1,iblk-1 - irwofs=irwofs+iblsz(i,2) - 1 continue - iclofs=iaclof - else - irwofs=iblsz(1,2) - iclofs=iblsz(1,1) - do 2 i=iscen+2,iblk-1 - irwofs=irwofs+iblsz(i,2) - iclofs=iclofs+iblsz(i,1) - 2 continue - endif -c - IF (FEASBL .NE. 'FEASIBLE') stop - CALL FREEVA(iblk,iclofs,irwofs,PROBLM,NFREE,VRFREE,ORGIDX,C) - CALL ISLACK(iblk,irwofs,PROBLM,VRFREE,ORGIDX,NFREE, - 1 C,ZEROD,NROWS) -C -C Since finding a P1 solution doesn't require modifying the A or AT data -C structures, I'm moving CSLACK and SETPH1 outside the loop that builds -C A,AT,AAT, etc. THen, the subs will be modified to be consistent with -C the new data structure. -C - 5 CONTINUE - IF(IBLK.GE.2.AND.IBLK.LE.ISCEN+1) THEN - IBLAAT(IBLK+1)=NIAAT - IDIAG(IBLK+1)=IDIAG(IBLK) - IBLPTR(IBLK+1)=IBLPTR(IBLK) - IBLLST(IBLK+1)=IBLLST(IBLK) - IBLDST(IBLK+1)=IBLDST(IBLK) - GOTO 7 - ENDIF - CALL STJAAT(iblk,ZEROI) -c do 195 i=1,ncols -c if(ja(i).le.0) write(6,*) ' ja(',i,')=',jaat(i) -c195 continue -c do 193 i=IBLAAT(iblk+1),IBLAAT(IBLK+2)-1 -c do 194 j=iaat(i),iaat(i+1)-1 -c if(jaat(j).le.0) write(6,*) ' jaat(',j,')=',jaat(j) -c194 continue -c193 continue - CALL OUTPRD(iblk,ERROR) - IF (ERROR .GT. ZEROI) then - write(0,*) 'after OUTPRD, error = ',error - stop - endif -C -C - CALL BLDMTR(iblk,iclofs,IRWOFS,D,ZEROD,EPSILO) -C - 7 CONTINUE -C -C - 15 CONTINUE -C All code previous deb'd -C ........................................ END FORVECION LOOP ............. -C -C - CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) -C - IF (SLKMIN .GT. ZEROD) GOTO 42 - PHASE=1 - CALL SETPH1(B,Y,SLKMIN,NROWS,BIGM,BIGMIU) - CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) -C /* If in P1, need to recalc D*/ -C /* Since BLDMTR used old slacks, update*/ - CALL UPDMTR(D,SLACK,-1.0D0) - WRITE(6,*) ' PHASE 1 NECESSARY, MIN SLACK=',SLKMIN - 42 CONTINUE -C ........................................ BEGIN SOLUTION LOOP ............ -C - STUPTM= DTIME(0) - itim(9)=0 -C - DO 100 ITERTN=1,100 -C - CALL FINDDY(ITERTN,B,DY) -C - itim(9)=itim(9)+itim(8)-itim(7) -C /* Following: impdrt,findv*/ - IF(PHASE.EQ.1) THEN - CALL IMPDRT(DY,V,W,D,B) - ENDIF -c write(30,*) ' ITERTN=',ITERTN -c write(30,909) (i,dy(i),i=1,nrows) -c909 format(5(' dy(',i4,')=',E13.7)) -C /* Now form Atdy */ -C /* must transfer dt to tmp*/ - CALL ATYCSL(2,DY,DV,PHASE,SLKMIN) -C - CALL STEPSZ(STEP,UBSTEP,GAMMA,SLACK,DV,NCOLS,STOPC) -C - if(itertn.eq.1) write(0,*) ' .............', - 1 ' Iteration Log ..............' - if(MOD(itertn,5).eq.0) then - if(STOPC.ge.99) write(0,1192) itertn, phase - if(STOPC.lt.99) write(0,1191) itertn, stopc, - 1 phase - endif - 1191 format(' Iteration: ',i5,' Stop Criterion: ',F15.12, - 1 ' Phase: ',I2) - 1192 format(' Iteration: ',i5,' Phase: ',i2) -C - IF ((STEP/GAMMA) .EQ. UBSTEP) THEN - FEASBL= 'STUCK-1 ' - GOTO 120 - END IF -C - DO 80 I= 1, NROWS - Y(I)= Y(I) + STEP*DY(I) - 80 CONTINUE - IF (PHASE .EQ. 2) GOTO 90 - Y(NROWS+1)= Y(NROWS+1) + STEP*DY(NROWS+1) - IF (Y(NROWS+1) .GT. EPSILO) GOTO 85 - IF (Y(NROWS+1) .GT. 0.0D0) GOTO 90 - ITERA1= ITERTN - PHS1TM= DTIME(0) -C - PHASE=2 -C - GOTO 90 -C - 85 IF ((YOLD - Y(NROWS+1)) .GE. EPSPH1) THEN - YOLD=Y(NROWS+1) - GOTO 95 - ENDIF - FEASBL= 'STUCK-2 ' -C /* Calc DUALOB */ - ZLOWER= 0.0D0 - DO 86 I=1,NROWS - ZLOWER= ZLOWER + B(I) * Y(I) - 86 CONTINUE - GOTO 120 -C /* Calc DUALOB */ - 90 ZLOWER= 0.0D0 - DO 91 I=1,NROWS - ZLOWER= ZLOWER + B(I) * Y(I) - 91 CONTINUE - IF (ZLOWER .GT. UBOUND) THEN - FEASBL= 'UNBOUNDD' - GOTO 120 - ENDIF -c -c A note regarding stopping criterion: Here, we use the criterion -c that the obj function does not change in the 8th sig. digit. -c To do comparisons w/ OSL, need to decrease RMUlimit to 0. -c Note: Max. comp. slackness doesn't work (alg. not numerically -c stable enough.) (i.e. max_{i} (x(i)*slack(i)) < stpeps -c - AUX=DABS(ZOLD) - IF (AUX .LT. 1.0D0) AUX=1.0D0 - STOPC=DABS(ZLOWER-ZOLD)/AUX - IF (STOPC .LE. EPSSTP) GOTO 120 - ZOLD=ZLOWER -C /* Find new slacks */ - 95 CALL ATYCSL(1,Y,ATY,PHASE,SLKMIN) -C - CALL UPDMTR(D,SLACK,EPSUPD) - - 100 CONTINUE -C - 120 IF (PHASE .EQ. 2) THEN - ITERA2= ITERTN-ITERA1 - PHS2TM= DTIME(0) - ELSE - ARTFCL= Y(NROWS+1) - ITERA1= ITERTN - PHS1TM= DTIME(0) - itim(5)=itim(3) - ENDIF -C /* Find Primal Objective */ - ZEROD =0.0D0 - ZUPPER=ZEROD - DO 122 I=1,NCOLS - X(I)= DV(I)/(SLACK(I)*SLACK(I)) - ZUPPER= ZUPPER+C(I)*X(I) - 122 CONTINUE -C /* GetZ value */ - IF (PROBLM .EQ. 'DUAL') GOTO 125 - DUMMY= ZUPPER - ZUPPER= -ZLOWER - ZLOWER= -DUMMY -C - 125 ZUPPER= ZUPPER + ZCONST - ZLOWER= ZLOWER + ZCONST -C - CALL GETSTX(X,NCOLS,ZEROI,NROWS,Y,PROBLM) -C 160 CALL PRINT(ZUPPER,ZLOWER,FEASBL,PROBLM,NFIX,NRWDRP, -C 1 NSTFRE,PHASE,ARTFCL,BIGM) -C - write(0,167) NZTOTL, ZLOWER - write(0,168) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), - 1 itim(5)-itim(4),itim(6)-itim(5),itim(9) - write(0,169) itim(6)-itim(1) - write(6,167) NZTOTL, ZLOWER - write(6,168) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), - 1 itim(5)-itim(4),itim(6)-itim(5),itim(9) - write(6,169) itim(6)-itim(1) -c write(6,170) itim(2)-itim(1),itim(3)-itim(2),itim(4)-itim(3), -c 1 itim(5)-itim(4),itim(6)-itim(5),itim(9),itim(6)-itim(1),NIAAT, - 2 ZLOWER - - 167 format(/' Nzs in AAT :',i15,' Solution :',E15.9) - 168 format(' Reading Time :',i15,' Formation Time :',i15/ - 1 ' Phase I Time :',i15,' Phase II Time :',i15/ - 2 ' Print Time :',i15,' TOTAL finddy time :',i15) - 169 format(' TOTAL TIME: : ',i15) - 170 format('SUM: ',6(I12,','),',,',2I12,E15.9) - 180 STOP - END ##### sed > BQ_Decomp/makefile <<'#####' 's/^-//' -.SUFFIXES : .o $(.SUFFIXES) - -PROBS = SC205 SCAGR7 SCFXM1 SCRS8 SCSD8 SCTAP - -OBJECTS = readin.o utils.o finddylin.o impdrt.o numerics.o wrapup.o blkio.o genmat.o main.o - -.f.o : - xlf -c -O $< - -bq2 : $(OBJECTS) timer.o - xlf -O -o bq2 $(OBJECTS) timer.o /usr/um/math/lib/linpack_d.a -L/afs/engin.umich.edu/group/engin/ioe/spk/AIX -lspk - -uncomp: - uncompress bq1.Z - cat bq1 | tar -xvf - - -all: bq2 ##### sed > BQ_Decomp/numerics.f <<'#####' 's/^-//' -C -C *************************************************************** -C * * -C * SUBROUTINE ATYCSL: This subroutine computes the slack * -C * in each constraint of the LP program * -C * for the current solution y, and * -C * returns the minimum value, OR computes* -C * Aty, according to the following * -C * * -C * * -C * IOPT=1: Find slack, returned through COMMON * -C * iOPT=2: Find At(y), returned through PARAMETER aty * -C * * -C * subroutines called: NONE -MULATY incorporated here * -C * * -C * called by: MAIN PROGRAM * -C * * -C *************************************************************** -C - SUBROUTINE ATYCSL(IOPT,Y,ATY,PHASE,SLKMIN) - INCLUDE (sizes.h) -C - DOUBLE PRECISION A,B,C,Y(SIZEVEC),SLACK, ATY(SIZEMAT) - DOUBLE PRECISION ZCONST,SLKMIN -C - INTEGER NROWS,NCOLS,IA,JA,PHASE, IOPT -C - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /INPT2/ B(SIZEVEC), C(SIZEMAT), ZCONST - COMMON /BODY/ SLACK(SIZEVEC) -C -C Form A^t.y -C /* NIA - NO OF COLS */ - IF(IOPT.LE.0.OR.IOPT.GE.3) THEN - WRITE(6,*) ' INTERNAL ERROR: (atycsl) IOPT=',IOPT - STOP - ENDIF -C - DO 5 I=1,NIA - ATY(I)=0.0D0 - 5 CONTINUE -C - DO 10 I=IBLA(1),IBLA(2)-1 - I1=I-IBLA(1)+1 - DO 9 J=IA(I),IA(I+1)-1 - ATY(I1)=ATY(I1)+Y(JA(J))*A(J) - 9 CONTINUE - 10 CONTINUE - J2=IBLSZ(1,2) - DO 20 I=2,ISCEN+1 - DO 19 J=IBLA(I),IBLA(I+1)-1 - J1=(J-IBLA(I)+1)+IACLOF - DO 18 K=IA(J),IA(J+1)-1 - ATY(J1)=ATY(J1)+A(K)*Y(J2+JA(K)) - 18 CONTINUE - 19 CONTINUE - J2=J2+IBLSZ(I,2) - 20 CONTINUE -C /* Now do Wty. I3=COL_OFS */ - I1=IBLSZ(1,2) - I2=IBLSZ(1,1) - DO 30 I=ISCEN+2,2*ISCEN+1 - DO 29 J=IBLA(I),IBLA(I+1)-1 - I3=I2+(J-IBLA(I)+1) - DO 28 K=IA(J),IA(J+1)-1 - ATY(I3)=ATY(I3)+A(K)*Y(I1+JA(K)) - 28 CONTINUE - 29 CONTINUE - I1=I1+IBLSZ(I,2) - I2=I2+IBLSZ(I,1) - 30 CONTINUE -C /* Now, i2= total no of cols*/ - IF(I2.NE.NCOLS) THEN - WRITE(6,*) ' IN CSLACK, I2,NCOLS=',I2,NCOLS - WRITE(6,*) ' STOPPING CONDITION ' - STOP - ENDIF -C /* Following was in MULATY */ - IF(PHASE.EQ.1) THEN - DO 35 I=1,I2 - ATY(I)=ATY(I)-Y(NROWS+1) - 35 CONTINUE - ENDIF -C /*If IOPT=2 then skip to end*/ - IF(IOPT.EQ.2) GOTO 50 - SLKMIN=1.0E+05 - DO 40 I=1,i2 - SLACK(I)= C(I) - ATY(I) - IF(SLACK(I).LT.SLKMIN) SLKMIN=SLACK(I) - 40 CONTINUE -C -C The next command avoids having a small slack in the Phase1 -C problem. This is why we are rounding the value from -1.0--0.0 -C to -1.0. -C - IF ((SLKMIN.GT.-1.0).AND.(SLKMIN.LE.0.0)) SLKMIN=-1.0D0 -C - 50 RETURN - END -C -C ************************************************************ -C * SUBROUTINE INTLZE, sete the values of the parameters * -C * and the initial values of some variables. * -C * * -C * called by: Main program subroutines called: None * -C ************************************************************ -C - SUBROUTINE INTLZE(EPSILO,EPSUPD,ERROR,FEASBL,GAMMA,NFIX,NFREE, - 1 NLBD,NRWDRP,NSTVAR,PHASE,PROBLM,UBOUND, - 1 UBSTEP,YOLD,ZCONST,ZOLD,NCHSGN,STOPC,BIGMIU,EPSPH1) - INCLUDE (sizes.h) -C - CHARACTER*3 BLANK1 - CHARACTER*4 PROBLM - CHARACTER*8 BLANK2, FEASBL -C - DOUBLE PRECISION ZEROD, UBOUND, ZCONST, GAMMA, - 1 EPSILO, EPSUPD, LST2, STOPC, - 1 UBSTEP, YOLD, ZOLD, BIGMIU, EPSPH1 -C - REAL STUPTM, PHS1TM, PHS2TM, MXITNS, PREPTM, READTM -C - INTEGER ERROR, ZEROI, ITERA1, ITERA2, NFIX, NFREE, NLBD, NRWDRP, - 1 NSTVAR, PHASE, NCHSGN, PTR2, DST2 -C - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /OUTPR2/ IBLPTR(MAXBLK),IBLLST(MAXBLK),IBLDST(MAXBLK), - 1 LST2(SIZEPRD),DST2(SIZEPRD), PTR2(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, - 1 STUPTM, PHS1TM, PHS2TM, MXITNS - COMMON /TMPAAT/ NZTOTL,NIAAT -C - BLANK1= ' ' - BLANK2= ' ' - BIGMIU= 1.0D5 - EPSILO= 1.0D-06 - EPSPH1= EPSILO - EPSUPD= 0.1 - ZEROI = 0 - ZEROD = 0.0D0 - ERROR = ZEROI - FEASBL= 'FEASIBLE' - GAMMA = 0.95D0 - STOPC=999 - ITERA1= ZEROI - ITERA2= ZEROI - MXITNS= 50 - NFIX = ZEROI - NFREE = ZEROI - NCHSGN= ZEROI - NLBD = ZEROI - NZTOTL= ZEROI - NRWDRP= ZEROI - - NSTVAR= ZEROI - PHASE = 2 - PROBLM= 'DUAL' - STUPTM = ZEROD - PHS1TM = ZEROD - PHS2TM = ZEROD - UBSTEP= 1.0E+20 - UBOUND= 1.0E+20 - YOLD = UBOUND - ZCONST= ZEROD - ZOLD = -1.0E+20 -C - ibla(1)=1 - iblat(1)=1 - iblaat(1)=1 - idiag(1)=1 - iscen = zeroi - narows= zeroi - nacols= zeroi - nwrows= zeroi - nwcols= zeroi - NIAAT = 1 -C - IBLPTR(1)=1 - IBLLST(1)=1 - IBLDST(1)=1 -C - RETURN - END - -C -C ********************************************************* -C * * -C * ISLACK: VERIFIED * -C * * -C * Subroutine ISLACK, inserts a slack variable into an * -C * inequation to modify it to an equation. * -C * * -C * called by: MAIN PROGRAM subrout. called: UPTARR * -C ********************************************************* -C - SUBROUTINE ISLACK(iblk,irwofs,PROBLM,VRFREE,ORGIDX, - 1 NFREE,C,ZEROD,NROWS) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE - CHARACTER*4 PROBLM - CHARACTER*8 NAMERW, NAMECL -C - DOUBLE PRECISION COEFSK, ZEROD, C(SIZEVEC) -C - INTEGER NROWS, I, J, VRFREE(SIZEVEC), NFREE, ORGIDX(SIZEVEC) -C - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /debug/ idbg - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof -C -C Only add slack on a,w'S, Problem should always be dual -C - IF (PROBLM .EQ. 'DUAL') GOTO 30 - COEFSK= -1.0D0 - DO 20 I=irwofs,nrows - DO 10 J=1,NFREE - IF (ORGIDX(I) .EQ. VRFREE(J)) GOTO 20 - 10 CONTINUE - CALL UPTARR(iblk,COEFSK,i,C,ZEROD) - 20 CONTINUE - GOTO 100 -C - 30 IF(IBLK.GE.2.AND.IBLK.LE.ISCEN+1) GOTO 100 - DO 70 I=1,IBLSZ(IBLK,2) - COEFSK= 1.0D0 - I1=IRWOFS+I - IF (RTYPE(I1+1) .EQ. ' E ') GOTO 70 - IF (RTYPE(I1+1) .EQ. ' N ') GOTO 80 - IF (RTYPE(I1+1) .EQ. ' G ') GOTO 50 - IF (RTYPE(I1+1) .EQ. ' L ') GOTO 60 - WRITE(6,40) I1+1, NAMERW(I1+1) - 40 FORMAT('ERROR. There is NO TYPE E,N,G or L for row ',I4/ - 1 'Name = ',A8) - GOTO 100 -C - 50 COEFSK= -1.0D0 - WRITE(6,*) 'ERROR, SHOULDNT BE g ROWS AT THIS STAGE' - 60 CALL UPTARR(iblk,COEFSK,i,C,ZEROD) -C - RTYPE(I1+1)= ' E ' -C - 70 CONTINUE - GOTO 100 -C - 80 WRITE(6,90) - 90 FORMAT('WARNING. Adding SLACK to the objective function?') -C - 100 RETURN - END -C - SUBROUTINE LUBKSB(A,N,NP,INDX,B) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION A(NP,NP),INDX(N),B(N) -C -C - II=0 - DO 12 I=1,N - LL=INDX(I) - SUM=B(LL) - B(LL)=B(I) - IF (II.NE.0)THEN - DO 11 J=II,I-1 - SUM=SUM-A(I,J)*B(J) -11 CONTINUE - ELSE IF (SUM.NE.0.0D0) THEN - II=I - ENDIF - B(I)=SUM -12 CONTINUE - DO 14 I=N,1,-1 - SUM=B(I) - IF(I.LT.N)THEN - DO 13 J=I+1,N - SUM=SUM-A(I,J)*B(J) -13 CONTINUE - ENDIF - B(I)=SUM/A(I,I) -14 CONTINUE - RETURN - END -C -C - SUBROUTINE LUDCMP(A,N,NP,INDX,D) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - PARAMETER (NMAX=300,TINY=1.0E-10) - DIMENSION A(NP,NP),INDX(N),VV(NMAX) - D=1.0D0 - DO 12 I=1,N - AAMAX=0.D0 - DO 11 J=1,N - IF (DABS(A(I,J)).GT.AAMAX) AAMAX=DABS(A(I,J)) -11 CONTINUE - IF (AAMAX.EQ.0.D0) THEN - WRITE(6,*) 'g1 OR G2 IS A Singular matrix.' - STOP - ENDIF - VV(I)=1.0D0/AAMAX -12 CONTINUE - DO 19 J=1,N - IF (J.GT.1) THEN - DO 14 I=1,J-1 - SUM=A(I,J) - IF (I.GT.1)THEN - DO 13 K=1,I-1 - SUM=SUM-A(I,K)*A(K,J) -13 CONTINUE - A(I,J)=SUM - ENDIF -14 CONTINUE - ENDIF - AAMAX=0.0D0 - DO 16 I=J,N - SUM=A(I,J) - IF (J.GT.1)THEN - DO 15 K=1,J-1 - SUM=SUM-A(I,K)*A(K,J) -15 CONTINUE - A(I,J)=SUM - ENDIF - DUM=VV(I)*DABS(SUM) - IF (DUM.GE.AAMAX) THEN - IMAX=I - AAMAX=DUM - ENDIF -16 CONTINUE - IF (J.NE.IMAX)THEN - DO 17 K=1,N - DUM=A(IMAX,K) - A(IMAX,K)=A(J,K) - A(J,K)=DUM -17 CONTINUE - D=-D - VV(IMAX)=VV(J) - ENDIF - INDX(J)=IMAX - IF(J.NE.N)THEN - IF(A(J,J).EQ.0.0D0)A(J,J)=TINY - DUM=1.0D0/A(J,J) - DO 18 I=J+1,N - A(I,J)=A(I,J)*DUM -18 CONTINUE - ENDIF -19 CONTINUE - IF(A(N,N).EQ.0.0D0)A(N,N)=TINY - RETURN - END -C -C ************************************************************* -C * * -C * Subroutine NORM, finds the euclidean norm of a given * -C * vector. * -C * * -C * called by: subroutine STARTY(JSTART) * -C * subroutines called: NONE. * -C ************************************************************* -C - SUBROUTINE NORM(X,NORMA,XSIZE) - INCLUDE (sizes.h) -C - DOUBLE PRECISION X(SIZEVEC), NORMA -C - INTEGER XSIZE, I -C - NORMA= 0.0D0 - DO 10 I=1,XSIZE - NORMA= NORMA+X(I)*X(I) - 10 CONTINUE - NORMA= DSQRT(NORMA) -C - RETURN - END -C -C ************************************************************* -C * * -C * Subroutine STARTY(JSTART), finds an initial value of Y given by * -C * (Norm(C)/Norm(At*b))*b. * -C * .....Performed before *any* preprocessing * -C * * -C * STARTY=0: use above: JSTART=1: read in new y * -C * called by: Main Program * -C * subroutines called: MULATY, NORM * -C ************************************************************* -C - SUBROUTINE STARTY(JSTART) - INCLUDE (sizes.h) -C - DOUBLE PRECISION ATB(SIZEVEC), AUX, NORMA,TA,TB,TC,TAUX - DOUBLE PRECISION ATY(SIZEVEC), Y,X,SLACK -C - INTEGER TNCOLS,TNROWS,TIA,TJA,TSORT,NSORT,TORIDX,j1 - INTEGER JSTART - LOGICAL TEST -C - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt -C - CALL NORM(TC,NORMA,TNCOLS) - AUX= NORMA -c - DO 20 J=1,tncols - ATB(J)= 0.0D0 - DO 10 I=tIA(J),tIA(J+1)-1 - ATB(J)= ATB(J)+tA(I)*tb(tJA(I)) - 10 CONTINUE - 20 CONTINUE -C - CALL NORM(ATB,NORMA,tNCOLS) - IF (NORMA .EQ. 0.0D0) NORMA=1.0D0 - AUX= AUX/NORMA -C - IF(JSTART.NE.0) write(6,*) '$$$$jstarT = ',JSTART - DO 30 I=1,tNROWS - IF(JSTART.EQ.0) THEN - Y(I)=AUX*tB(I) - ELSE IF (JSTART.EQ.1) THEN - Y(I)=AUX - ELSE IF(JSTART.eq.2) THEN - Y(I)=1.0D0 - ELSE IF(JSTART.eq.3) THEN - Y(I)=0.0D0 - ELSE - WRITE(6,*) ' ERROR: IN STARTY, JSTART = ',JSTART - ENDIF -C - 30 CONTINUE -C /* fIND sLACK */ -C It's necessary to find the slacks here and then augment them with y's -C from slacks added to iblk=1 (A). The slacks will eventually be -C overwritten with the correct slacks as ISLACK is performed.... -C -C Following is less than efficient: change later -C only need to calc slack for (ATTTT)t portion - rest calc'd later -C /* First, not in Tsort */ - J2=0 - DO 40 J=1,nacols - TEST=.TRUE. - DO 33 J1=1,NSORT - IF(TSORT(J1).EQ.J) TEST=.FALSE. - 33 CONTINUE - IF(.not.TEST) GOTO 40 - ATY(J)= 0.0D0 - DO 35 I=tIA(J),tIA(J+1)-1 - ATY(J)= ATY(J)+tA(I)*Y(tJA(I)) - 35 CONTINUE - J2=J2+1 - SLACK(J2)=TC(J)-ATY(J) - if(abs(slack(j2)).le.1.D-02) slack(j2)=-1.D-02 - 40 CONTINUE -C - DO 50 J=1,nacols -c IF(J.GT.IMXSRT) GOTO 50 - TEST=.TRUE. - DO 43 J1=1,NSORT - IF(TSORT(J1).EQ.J) TEST=.FALSE. - 43 CONTINUE - IF(TEST) GOTO 50 - ATY(J)= 0.0D0 - DO 45 I=tIA(J),tIA(J+1)-1 - ATY(J)= ATY(J)+tA(I)*Y(tJA(I)) - 45 CONTINUE - J2=J2+1 - SLACK(J2)=TC(J)-ATY(J) - if(abs(slack(j2)).le.1.D-02) slack(j2)=-1.D-02 - 50 CONTINUE -C - RETURN - END -C -C -C ************************************************************ -C * * -C * Subroutine STEPSZ : calculates the step size. The * -C * next point Y is obtained by * -C * moving a step size from the current point * -C * in the direction dy. * -C * * -C * called by: MAIN program subroutines called: NONE * -C ************************************************************ -C - SUBROUTINE STEPSZ(STEP,UBSTEP,GAMMA,SLACK,DV,NCOLS,STPCRI) - INCLUDE (sizes.h) -C - DOUBLE PRECISION STEP,GAMMA,UBSTEP,SLACK(SIZEVEC) - DOUBLE PRECISION STPCRI, TEMP, DV(SIZEMAT) -C - INTEGER I,NCOLS -C - STEP= UBSTEP -c STPCRI=0.0D0 - DO 10 I=1,NCOLS - IF (DV(I) .LE. 0.0D0) GOTO 10 - TEMP= SLACK(I)/DV(I) - IF (TEMP .LT. STEP) STEP=TEMP - 10 CONTINUE - STEP= GAMMA*STEP -C - 40 RETURN - END ##### sed > BQ_Decomp/readin.f <<'#####' 's/^-//' - -C *********************************************************** -C * * -C * Program READAT: reads in the data of a LP program * -C * presented in standard MPS format, set the * -C * problem in the standard form and feeds in * -C * the appropriate arrays of the affine * -C * scaling program (Karmarkar's algorithm) * -C * to solve this LP. * -C * * -C * Subroutines called: READRW, READCL, READRH, READRG, * -C * READBD, ISLACK * -C * * -C * Written by: Jose' Arantes Jan/89 * -C *********************************************************** -C - SUBROUTINE READAT(JSTART,EPSPH1,EPSILO,BIGMIU) - INCLUDE (sizes.h) -C - CHARACTER*3 NEXTLN, tRTYPE, BLANK1 - CHARACTER*8 tNAMER, BLANK2, tNAMEC - character*20 argv -C - DOUBLE PRECISION tA, tB, tC, VALFIX, VALLBD, ZEROD, taux - DOUBLE PRECISION EPSILO, BIGMIU, EPSPH1 - INTEGER I, tNROWS, tNCOLS, ZEROI, NRWDRP, NFIX, NLBD, - 1 ORGIDX, VARFIX, VARLBD, NSTVAR, VRFREE, NFREE, - 1 VRCSGN, NCHSGN,tia, tja, tsort, nsort, TORIDX, - 1 JSTART -C - COMMON /tCONST/ tNAMER(SIZEVEC),tRTYPE(SIZEVEC),tNAMEC(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /BLANK/ ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /FIXED/ VALFIX(SIZEBLK),VALLBD(SIZEBLK),VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK), VRFREE(SIZEVEC), VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - I= ZEROI - tncols=zeroi - tnrows=zeroi - JSTART=0 - 10 NEXTLN= BLANK1 - READ(4,15) NEXTLN - 15 FORMAT(A3) -C - IF (NEXTLN .EQ. 'ROW') GOTO 20 -C - call getarg(2,argv) - read(argv,*) iscen - call getarg(3,argv) - read(argv,*) nacols - call getarg(4,argv) - read(argv,*) narows - call getarg(5,argv) - read(argv,*) nwcols - call getarg(6,argv) - read(argv,*) nwrows - if(nextln.eq.'STO') then - read(4,*) iscen, nacols, narows, nwcols, nwrows - endif - if(nextln.eq.'PAR') then - read(4,*) bigmiu, epsilo,epsph1 - else - bigmiu=1000000 - epsilo=.000001 - epsph1=.000001 - endif - if(nextln.eq.'STA') JSTART=1 -C - IF (I .EQ. 50) GOTO 30 - GOTO 10 -C - 20 I= ZEROI - CALL READRW(NEXTLN) - IF (NEXTLN .NE. 'COL') GOTO 50 - CALL READCL(NEXTLN) - IF (NEXTLN .NE.'RHS') GOTO 70 - CALL READRH(NEXTLN) -c IF (NEXTLN .EQ. 'RAN') CALL READRG(NEXTLN) -c IF (NEXTLN .EQ. 'BOU') CALL READBD(NEXTLN) - IF (NEXTLN .EQ. 'END') GOTO 110 - GOTO 90 -C - 30 WRITE(6,40) - 40 FORMAT('ERROR: The indicator card ROWS was not found.') - GOTO 110 - 50 WRITE(6,60) - 60 FORMAT('ERROR: The indicator card COLUMNS was not found.') - GOTO 110 - 70 WRITE(6,80) - 80 FORMAT('ERROR: The indicator card RHS was not found.') - GOTO 110 - 90 WRITE(6,100) - 100 FORMAT('ERROR: The indicator card ENDDATA was not found.') -C - 110 CONTINUE - RETURN - END -C -C ********************************************************** -C * * -C * Subroutine READCL: it reads in the cost coefficients * -C * and the non-zero elements of each column, associated * -C * row index and value. It also fills in the arrays * -C * C(I), IA(I), JA(I) and A(I), where: * -C * * -C * C(I) has the cost coefficient of the ith variable * -C * IA(I) has the value of the position in JA and A * -C * where starts the Ith column. * -C * JA(K) if the current column is J, then JA(K) gives * -C * the row index of the (K-IA(J)+1)th non-zero * -C * element in column J. * -C * A(K) similarly A(K) has the value of the (K-IA(J)+1) * -C * th non-zero element in column J. * -C * * -C * called by: program READAT * -C * Subroutines called : INDROW * -C ********************************************************** -C - SUBROUTINE READCL(NEXTLN) - INCLUDE (sizes.h) -C - CHARACTER*3 NEXTLN, BLANK1, tRTYPE - CHARACTER*8 CURRCL, LASTCL, ROW1, ROW2, BLANK2, tNAMEC - CHARACTER*8 tNAMER -C - DOUBLE PRECISION A1, A2, tA, tc, tB, ZEROD, VALFIX, - 1 VALLBD, TAUX -C - INTEGER tNCOLS, tNROWS, I, INDR, ZEROI, NFIX, NLBD, - 1 NSTVAR, ORGIDX, VRFREE, NFREE, NRWDRP, VARFIX, - 1 VRCSGN, NCHSGN, VARLBD, tia, tja, TSORT, NSORT, - 1 toridx -C - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) - COMMON /BLANK/ ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /FIXED/ VALFIX(SIZEBLK), VALLBD(SIZEBLK), VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK), VRFREE(SIZEVEC), VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - DO 5 I=1,3000 - tC(I)= ZEROD - tNAMEC(I)= BLANK2 - 5 CONTINUE - LASTCL= 'STARTCOL' - tNCOLS= ZEROI - I= ZEROI -C - 10 CURRCL= BLANK2 - NEXTLN= BLANK1 - ROW1= BLANK2 - ROW2= BLANK2 - A1= ZEROD - A2= ZEROD -C - READ(4,20) NEXTLN,CURRCL,ROW1,A1,ROW2,A2 - 20 FORMAT(A3,1X,A8,2X,A8,2X,F12.4,3X,A8,2X,F12.4) -C - IF (NEXTLN .EQ. 'RHS') GOTO 70 - IF (CURRCL .EQ. LASTCL) GOTO 30 - tNCOLS= tNCOLS+1 - TORiDX(tNCOLS)= tNCOLS - tIA(tNCOLS)= I+1 - LASTCL= CURRCL - tNAMEC(tNCOLS)= CURRCL -C - 30 CALL INDROW(INDR,ROW1,tNROWS) -C The program ignores coefic. if row is not defined (INDR=9999) - IF (INDR .EQ. 9999) GOTO 50 - IF (INDR .NE. ZEROI) GOTO 40 - tc(tNCOLS)= A1 - GOTO 50 - 40 I= I+1 - tJA(I)= INDR -C -C note that the correct sign is applied here: don't have to worry -C about it later -C - tA(I)= A1*taux(indr) - 50 IF (ROW2 .EQ. BLANK2) GOTO 10 - CALL INDROW(INDR,ROW2,tNROWS) - IF (INDR .EQ. 9999) GOTO 10 - IF (INDR .NE. ZEROI) GOTO 60 - tc(tNCOLS)= A2 - GOTO 10 - 60 I=I+1 - tJA(I)= INDR - tA(I)= A2*taux(indr) - GOTO 10 -C ------------------------------------------------- -C IA(NCOLS+1) points to where finishes the non-zero entries of A -C -------------------------------------------------------------- - 70 tIA(tNCOLS+1)= I+1 - NSTVAR= tNCOLS -C - write(6,*) ' No. of Non-Zeros:',i - write(6,*) ' Density :',float(i)/float(tnrows*tncols) - write(6,*) ' NCOLS, NROWS :',tncols, tnrows - RETURN - END -C -C ********************************************************** -C * * -C * Subroutine READRH: it reads in the RHS associated * -C * to each constraint of the LP program. * -C * * -C * Assumption: Only one RHS vector is provided in * -C * the MPS format. * -C * * -C * called by: Program READAT * -C * Subroutines called: INDROW * -C ********************************************************** -C - SUBROUTINE READRH(NEXTLN) - INCLUDE (sizes.h) -C - CHARACTER*3 NEXTLN, BLANK1 - CHARACTER*8 ROW1,ROW2,CURRHS, BLANK2 -C - DOUBLE PRECISION A1, A2, ZEROD, tb, tc, taux, ta -C - INTEGER tNROWS, I, INDR, ZEROI, TSORT, NSORT,tncols, tia, tja, - 1 toridx -C - COMMON /BLANK/ ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt -C - DO 10 I=1,tNROWS - tB(I)= ZEROD - 10 CONTINUE -C - 20 ROW1= BLANK2 - ROW2= BLANK2 - A1= ZEROD - A2= ZEROD - NEXTLN= BLANK1 -C - READ(4,30) NEXTLN, CURRHS, ROW1, A1, ROW2, A2 - 30 FORMAT(A3,1X,A8,2X,A8,2X,F12.4,3X,A8,2X,F12.4) -C - IF (NEXTLN .NE. BLANK1) GOTO 40 - CALL INDROW(INDR,ROW1,tNROWS) - tB(INDR)= A1*TAUX(INDR) - IF (ROW2 .EQ. BLANK1) GOTO 20 - CALL INDROW(INDR,ROW2,tNROWS) - tB(INDR)= A2*TAUX(INDR) - GOTO 20 -C - 40 RETURN - END -C -C ********************************************************** -C * * -C * Subroutine READRW: it reads in the row name and the * -C * (in) equality sign associated to each row. * -C * The signs G, L,E correspond to constraints * -C * and N to the objective function. * -C * * -C * called by: Program READAT subroutines called: None * -C * * -C ********************************************************** -C - SUBROUTINE READRW(NEXTLN) - INCLUDE (sizes.h) -C - CHARACTER*3 NEXTLN,tRTYPE, BLANK1, AUX1 - CHARACTER*8 tNAMER, BLANK2, tNAMEC, AUX2 -C - DOUBLE PRECISION ZEROD, TB, TC,TAUX,ta -C - INTEGER I, J, ZEROI,tnrows,TSORT,NSORT,tia,tja,toridx,tncols - LOGICAL TEST -C - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) - COMMON /BLANK/ ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt -C - I= 1 - 10 tRTYPE(I)= BLANK1 - tNAMER(I)= BLANK2 - READ(4,20) tRTYPE(I), tNAMER(I) -C WRITE(6,*) ' READING ROW: ',TNAMER(I) - 20 FORMAT(A3,1X,A8) - IF (tRTYPE(I) .EQ. 'COL') GOTO 30 - I= I+1 - GOTO 10 -C - 30 tNROWS= I-1 - NEXTLN= 'COL' - IF (tRTYPE(1) .EQ. ' N ') TEST=.TRUE. - DO 40 J=2,tNROWS - tAUX(j-1)=1.0D0 - IF (tRTYPE(j) .eq. ' G ') THEN - tRTYPE(j)= ' L ' - tAUX(j-1)= -1.0D0 - END IF - IF (tRTYPE(J) .EQ. ' N '.AND..NOT.TEST) THEN - AUX1= tRTYPE(J) - AUX2= tNAMER(J) - tNAMER(J)= tNAMER(1) - tRTYPE(J)= tRTYPE(1) - tNAMER(1)= AUX2 - tRTYPE(1)= AUX1 - TEST=.TRUE. - ENDIF - 40 CONTINUE - IF(.NOT.TEST) THEN - WRITE(6,*) 'NO OBJECTIVE FUNCTION !!!!' - STOP - ENDIF -C ---------------------------------- -C The obj. funct. will always be the 1st row in NAMERW,RTYPE -C ---------------------------------------------------------- -C - 60 tNROWS= tNROWS-1 -C - 70 RETURN - END ##### sed > BQ_Decomp/sizes.h <<'#####' 's/^-//' - INTEGER SIZEBLK,SIZEMAT,SIZEAAT,SIZEPRD,SIZEDSAVE, - 1 MAXBLK,SIZEVEC,SIZET - DOUBLE PRECISION EPSSTP - PARAMETER (SIZEBLK = 14000) - PARAMETER (SIZET = 300) - PARAMETER (SIZEVEC = 12000) - PARAMETER (SIZEMAT = 40000) - PARAMETER (SIZEAAT = 160000) - PARAMETER (SIZEPRD = 160000) - PARAMETER (SIZEDSAVE = 200000) - PARAMETER (MAXBLK = 520) - PARAMETER (EPSSTP = 5.0D-08) ##### sed > BQ_Decomp/timer.c <<'#####' 's/^-//' -#include -#include -#include -/* Options: 1 = intitialize, 2= get time, difftime*/ - -int dfhtime(int *iopt) -{ - static time_t t1, t2; - if(*iopt==1) { - t1=time(NULL); - printf(" TIME: Initialize: %s",ctime(&t1)); - } - else { - t2=time(NULL); - printf(" TIME: %s :",ctime(&t2)); - fprintf(stderr," DIFF: %f secs.\n",difftime(t2,t1)); - printf(" DIFF: %f secs.\n",difftime(t2,t1)); - t1=t2 ; - } - return 0; -} - ##### sed > BQ_Decomp/utils.f <<'#####' 's/^-//' -C -C *********************************************************** -C * Subroutine COMPCT: sorts IAUX in crescent order and * -C * eliminates duplicated entries. * -C * * -C * called by: STJAAT subroutines called: NONE * -C *********************************************************** -C - SUBROUTINE COMPCT(IAUX,TOTAL) - INCLUDE (sizes.h) -C - INTEGER IAUX(400), TOTAL, TTSEEN, IBEG, IEND, IDXMIN, MIN, I -C - TTSEEN= 1 - 10 IF (TTSEEN .GE. TOTAL) GOTO 40 - MIN= IAUX(TTSEEN) - IDXMIN= TTSEEN - IBEG= TTSEEN+1 - IEND= TOTAL -C - DO 30 I= IBEG, IEND - 20 IF (MIN .LT. IAUX(I)) GOTO 30 - IF (MIN .EQ. IAUX(I)) THEN - IAUX(I)= IAUX(TOTAL) - IAUX(TOTAL)= 99999 - TOTAL= TOTAL-1 - GOTO 20 - ENDIF - MIN= IAUX(I) - IDXMIN= I - 30 CONTINUE -C - IAUX(IDXMIN)= IAUX(TTSEEN) - IAUX(TTSEEN)= MIN - TTSEEN= TTSEEN+1 - GOTO 10 -C - 40 RETURN - END -C -C ********************************************************* -C * * -C * DROPCL NOT YET VERIFIED * -C * * -C * Subroutine DROPCL, eliminates a column of the LP by * -C * updating NCOLS and the arrays C, NAMECL, * -C * IA, JA, A AND ORGIDX. * -C * * -C * called by: SUBROUTINES FIXVAR and FREEVA * -C * subroutines called: None. * -C ********************************************************* -C - SUBROUTINE DROPCL(iblk,COLUMN) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE, BLANK1 - CHARACTER*8 NAMERW, NAMECL, BLANK2 -C - DOUBLE PRECISION A, B, C, VALFIX, VALLBD, ZCONST, ZEROD -C - INTEGER COLUMN, NROWS, NCOLS, IA, JA, ZEROI, IBEG, IEND, - 1 I, TOTNZC, VARFIX, VARLBD, NFIX, NLBD, NSTVAR, - 1 ORGIDX, VRFREE, NFREE, NRWDRP, VRCSGN, NCHSGN, - 1 ofstcl, itotc,itotr,nja,nia,niat,njat,iaclof -C -C ofstcl=offsetted column -C - COMMON /debug/ idbg - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /INPT2/ B(SIZEVEC), C(SIZEMAT), ZCONST - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /FIXED/ VALFIX(SIZEBLK), VALLBD(SIZEBLK), VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK), VRFREE(SIZEVEC), VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - ofstcl=ibla(iblk)+column-1 - TOTNZC= IA(ofstcl+1)-IA(ofstcl) - NCOLS= NCOLS-1 -C - IF (ofstcl .EQ. (NCOLS+1)) GOTO 25 - DO 10 I=ofstcl,Nia - C(I)= C(I+1) - NAMECL(I)= NAMECL(I+1) - ORGIDX(I)= ORGIDX(I+1) - IA(I+1)= IA(I+2) - TOTNZC - 10 CONTINUE -C - ibla(iblk+1)=ibla(iblk+1)-1 -C - IBEG= IA(ofstcl) - IEND= IA(NCOLS+1)-1 - DO 20 I=IBEG,IEND - JA(I)= JA(I+TOTNZC) - A(I)= A(I+TOTNZC) - 20 CONTINUE -C - 25 IBEG= IA(NCOLS+1) - IEND= IA(NCOLS+2) - DO 30 I= IBEG,IEND - JA(I)= ZEROI - A(I)= ZEROD - 30 CONTINUE -C - IA(NCOLS+2)= ZEROI - C(NCOLS+1)= ZEROD - NAMECL(NCOLS+1)=BLANK2 -C - nja=nja-totnzc - nia=nia-1 - - RETURN - END -C -C ********************************************************* -C * * -C * * -C * FIXVAR NOT YET VERIFIED * -C * * -C * * -C * Subroutine FIXVAR: it eliminates a column of the LP * -C * which corresponds to a variable whose * -C * value is fixed. This implies, update RHS, * -C * update ZCONST and drop a column. * -C * * -C * called by: subroutine READBD, PREPRO, ONENZR * -C * subroutines called: INDCOL, DROPCL * -C ********************************************************* -C - SUBROUTINE FIXVAR(iclofs,irwofs,VALUE,CLNAME,iblk) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE,BLANK1 - CHARACTER*8 CLNAME, NAMERW, NAMECL, BLANK2 -C - DOUBLE PRECISION VALUE, A, B, C, VALFIX, VALLBD, ZCONST, - 1 ZEROD -C - INTEGER IDX, IA, JA, NROWS, NCOLS, I, IBEG, IEND, ZEROI, - 1 NFIX, NFREE, NLBD, NSTVAR, VARFIX, VARLBD, ORGIDX, - 1 VRFREE, NRWDRP, NCHSGN, VRCSGN, ibla, iblat, iblaat, - 1 idiag -C - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /INPT2/ B(SIZEVEC), C(SIZEMAT), ZCONST - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /FIXED/ VALFIX(SIZEBLK), VALLBD(SIZEBLK), VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK), VRFREE(SIZEVEC), VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - CALL INDCOL(IDX,CLNAME,iblk,iclofs) - NFIX= NFIX+1 - VARFIX(NFIX)= ORGIDX(IDX) - VALFIX(NFIX)= VALUE -C ---------------- -C UPDATE RHS -C ------------------ - IBEG= IA(ibla(iblk)+IDX-1) - IEND= IA(ibla(iblk)+IDX)-1 - DO 10 I=IBEG,IEND - B(JA(I)+irwofs)= B(JA(I)+irwofs) - A(I)*VALUE - 10 CONTINUE -C --------------------------------------------- -C Update the constant of the objective function -C ---------------------------------------------- - ZCONST= ZCONST + C(IDX+iclofs)*VALUE -C - CALL DROPCL(iblk,IDX) -C - RETURN - END -C -C -C ********************************************************** -C * * -C * FREEVA NOT YET VERIFIED * -C * * -C * Subroutine FREEVA, finds free variables and replace * -C * them by the difference of two non-negative * -C * variables. * -C * * -C * called by: Main Program subroutines called:FREVA1 * -C ********************************************************** -C - SUBROUTINE FREEVA(iblk,iclofs,irwofs,PROBLM,NFREE, - 1 VRFREE,ORGIDX,C) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE - CHARACTER*4 PROBLM - CHARACTER*8 NAMERW, NAMECL -C - DOUBLE PRECISION A, C(SIZEVEC), AT -C - INTEGER I, NCOLS, NFREE, VRFREE(SIZEVEC), ORGIDX(SIZEVEC), idiag, - 1 NROWS,IA, JA, J, iblk, ibla, iclofs, irwofs, iaclof, - 1 nia,nja,itotc,niat,njat, itotr, IAT, JAT, iblat, iblaat -C - COMMON /debug/ idbg - COMMON /CONSTR/ NAMERW(SIZEVEC),RTYPE(SIZEVEC),NAMECL(SIZEVEC) - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) -C -C PRoblem should be dual -C - IF (PROBLM .EQ. 'DUAL') GOTO 20 -C -C Using NCOLS should be ok, remember that we;ve only processed -C up to a point -C - DO 10 I=iclofs,ncols - IF (RTYPE(I+1) .NE. ' E ') GOTO 10 - CALL FREVA1(iblk,iclofs,irwofs,IA,JA,A,C,ncols, - 1 ibla(iblk)+(i-iclofs+1)) - 10 CONTINUE - GOTO 70 -C - 20 DO 60 I=1,NFREE - WRITE(6,*) 'WARNING, IN FREEVA, NFREE=',NFREE -C - DO 30 J=iclofs,ncols - if(idbg.ge.2) then - write(7,41) i,j,vrfree(i),orgidx(j) - endif - IF (VRFREE(I) .EQ. ORGIDX(J)) GOTO 50 - 30 CONTINUE - WRITE(6,40) VRFREE(I) - 40 FORMAT('The free variable',I4,' was not found in ORGIDX(J)') - 41 format(' FREEVA: i,j,vrfree(i),orgidx(j)',4(i4,2x)) - GOTO 70 -C - 50 CALL FREVA1(iblk,iclofs,irwofs,IA,JA,A,C,NCOLS, - 1 ibla(iblk)+(j-iclofs+1)) - ORGIDX(NCOLS)= ORGIDX(J) - 60 CONTINUE -C - 70 RETURN - END -C -C -C ********************************************************** -C * * -C * FREEVA1 NOT YET VERIFIED * -C * * -C * Subroutine FREVA1, creates a new column when replacing* -C * a free variable by the difference of two * -C * non-negative variables. Updates the arrays * -C * C, IA, JA, A, and the scalar NCOLS. * -C * * -C * Also, updates iat, jat, at ??????? * -C * * -C * HERE, JSTAR is big picture: (NOT 1..N?cols * -C * * -C * called by: subroutine FREEVA * -C * subroutines called: INDCOL * -C ********************************************************** -C - SUBROUTINE FREVA1(iblk,iclofs,irwofs,IA,JA,A,C,NCOLS,JSTAR) - INCLUDE (sizes.h) -C - DOUBLE PRECISION A(SIZEMAT), C(SIZEVEC) -C - INTEGER NCOLS, IA(SIZEVEC), JA(SIZEMAT), BEGCOL, ENDCOL, JSTAR, - 1ibla, iblat, iblaat,idiag, iblsz, itotr, nja,nia, - 1itotc,niat,njat,iaclof -C - COMMON /debug/ idbg - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof -C - if(idbg.ge.3) then - write(7,90) iblk,iclofs,irwofs,NCOLS,JSTAR - 90 format('FREVA1:iblk,iclofs,irwofs,NCOLS,JSTAR',4(2x,i5)) - endif - BEGCOL= IA(JSTAR) - ENDCOL= IA(JSTAR+1)-1 -C - IA(NCOLS+2)= IA(NCOLS+1)+ENDCOL-BEGCOL+1 -C - DO 10 I= BEGCOL,ENDCOL - nja=nja+1 - JA(nja)= JA(I) - A(nja)= -A(I) -C - 10 CONTINUE -C - NCOLS=NCOLS+1 - nia=nia+1 - iblsz(iblk,1)=iblsz(iblk,1)+1 - ibla(iblk+1)=ibla(iblk+1)+1 -C Update nacols - C(NCOLS)= -C(JSTAR) -C -C May have to update IAT, JAT, AT - WRITE(7,*) 'CALLED FREEVA1: SHOULD NEVER HAPPEN!!' -C - RETURN - END -C -C ********************************************************* -C * * -C * Subroutine DROPRW, it eliminates a constraint (a row)* -C * of the current LP by updating NROWS and * -C * the arrays RTYPE,NAMERW, B, JA, A, and IA * -C * * -C * called by: MAIN PROGRAM and subroutine PREPRO * -C * subroutines called: None. * -C ********************************************************* -C -C ********************************************************** -C * * -C * Subroutine INSERW: inserts a new constriant in the * -C * current LP. That is, adds one name in * -C * NAMERW(I) and one sign in RTYPE(I), and * -C * updates IA(I), JA(I), A(I) and B(I). * -C * * -C * called by: subroutines READRG, LOUPVA, SETPH1 * -C * subroutines called: None * -C ********************************************************** -C -C ********************************************************* -C * * -C * Subroutine CHSIGN: it changes the variable x with * -C * no lower bound and finite upper bound * -C * to y= -x with a finite lower bound and * -C * no upper bound. * -C * * -C * called by: subroutines READBD, PREPRO * -C * subroutines called: INDCOL. * -C ********************************************************* -C -C -C ********************************************************* -C * * -C * Subroutine INDCOL: it associates an index number * -C * 1, 2, ..n to each column name. * -C * * -C * called by: subroutine LOUPVA, LOBOVA, FIXVAR * -C * FREEVA, CHSIGN * -C * subroutines called: None. * -C ********************************************************* -C - SUBROUTINE INDCOL(INDEX,COLX,iblk,iclofs) - INCLUDE (sizes.h) -C - CHARACTER*3 tRTYPE - CHARACTER*8 tNAMEC, COLX,tNAMER -C - INTEGER INDEX, iblk, I, ncol -C - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) -C - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - ncol=nwcols - if(iblk.le.iscen+1) ncol=nacols - DO 10 I=1,NCOL - IF (tNAMEC(I+iclofs) .EQ. COLX) GOTO 30 - 10 CONTINUE -C - WRITE(6,20) COLX - 20 FORMAT(20X,'WARNING. Column index OF ',A8,' was not found.') - DO 25 I=1,70 - WRITE(6,24) tNAMEC(I) - 24 FORMAT('tNAMEC=',A8) - 25 CONTINUE - INDEX= 9999 - GOTO 40 -C - 30 INDEX= I -C - 40 RETURN - END - -C ********************************************************** -C * * -C * Subroutine INDROW: it associates an index number, * -C * 1,2 ..n to each constraint name. * -C * * -C * called by: subroutines READCL, READRH, READRG * -C * subroutines called: None. * -C ********************************************************** -C - SUBROUTINE INDROW(INDEX,ROWX,NROWS) - INCLUDE (sizes.h) -C - CHARACTER*3 tRTYPE - CHARACTER*8 tNAMER, ROWX, tNAMEC -C - INTEGER INDEX, NROWS, I -C - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) -C - DO 10 I=1,NROWS+1 - IF (tNAMER(I) .EQ. ROWX) GOTO 30 - 10 CONTINUE -C - WRITE(6,20) ROWX - 20 FORMAT(10X,'ROW INDEX NOT FOUND: TRYING TO FIND ',A8) - INDEX= 9999 - goto 40 -C - 30 INDEX= I-1 -C - 40 RETURN - END -C -C ************************************************************ -C * SUBROUTINE SORTCL, it sorts, inside each column, the * -C * nonzero coeficients in crescent * -C * order of row index "JA(I)". * -C * * -C * called by: Main Program subroutines called: None * -C ************************************************************ -C - SUBROUTINE SORTCL - INCLUDE (sizes.h) -C - DOUBLE PRECISION AUX, tA, tb, tc, taux -C - INTEGER tIA, tJA, BEGCL1, BEGCL2, ENDCL1, ENDCL2, I1, I2, MIN - INTEGER tNCOLS, tNROWS, J, NOZERO, IAUX, tsort - INTEGER TORIDX,narows,nwrows,nacols,nwcols -C - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof -C - icntr=0 - imxsrt=-999 - DO 30 J=1,tNCOLS - NOZERO= tIA(J+1)-tIA(J) - BEGCL1=tIA(J) - ENDCL1=tIA(J+1)-2 - ENDCL2=tIA(J+1)-1 - IF (NOZERO .LT. 2) GOTO 25 - DO 20 I1=BEGCL1,ENDCL1 - MIN=I1 - BEGCL2=I1+1 - DO 10 I2=BEGCL2,ENDCL2 - IF (tJA(I2) .GE. tJA(MIN)) GOTO 10 - MIN=I2 - 10 CONTINUE - AUX= tA(MIN) - IAUX= tJA(MIN) - tA(MIN)=tA(I1) - tJA(MIN)=tJA(I1) - tA(I1)= AUX - tJA(I1)=IAUX - 20 CONTINUE -C -C If beginnning of column is in A and End is in T, put in TSORT -C - 25 if(((tja(begcl1).le.narows).and.(tja(endcl2).gt.narows)) - 1 .or.((tja(begcl1).gt.narows).and.j.le.nacols)) then - icntr=icntr+1 - tsort(icntr)=j - if(imxsrt.le.j) imxsrt=j - endif - 30 CONTINUE -C - NSORT=icntr -C - 110 iblsz(1,1)=nacols - iblsz(1,2)=narows - do 120 i=1,iscen - iblsz(i+1,1)=nsort - iblsz(i+1,2)=nwrows - iblsz(iscen+i+1,1)=nwcols - iblsz(iscen+i+1,2)=nwrows - 120 continue - IACLOF=NACOLS-NSORT - RETURN - END -C -C Following are headers for unused subroutines for Birgeand Qi implementation -C For now, I'll assume that the input matrix doesn't need much preprocessing -C in the sense that there are no bounds/ranges, no vars are free, no rows, etc. -C are easily verified to be redundant. SO the only thing necessary is to get -C the problem in Ax=b,x>=0 form. (These are done by ISLACK). -C -C **************************************************************** -C * Subroutine PREPRO: Preprocessing the Input Data. Basically, * -C * * -C * eliminated * -C * * -C **************************************************************** -C -C ********************************************************* -C * * -C * Subroutine DROPRW, it eliminates a constraint (a row)* -C * of the current LP by updating NROWS and * -C * the arrays RTYPE,NAMERW, B, JA, A, and IA * -C * * -C * called by: MAIN PROGRAM and subroutine PREPRO * -C * subroutines called: None. * -C ********************************************************* -C -C ********************************************************** -C * * -C * Subroutine INSERW: inserts a new constriant in the * -C * current LP. That is, adds one name in * -C * NAMERW(I) and one sign in RTYPE(I), and * -C * updates IA(I), JA(I), A(I) and B(I). * -C * * -C * called by: subroutines READRG, LOUPVA, SETPH1 * -C * subroutines called: None * -C ********************************************************** -C -C ********************************************************* -C * * -C * Subroutine CHSIGN: it changes the variable x with * -C * no lower bound and finite upper bound * -C * to y= -x with a finite lower bound and * -C * no upper bound. * -C * * -C * called by: subroutines READBD, PREPRO * -C * subroutines called: INDCOL. * -C ********************************************************* -C ARWISE in another subroutine -C -C *********************************************************** -C * Subroutine ARWISE, stores the coeficient matrix A in * -C * a 'rowise' form through the arrays JAT, IAT, AT. * -C * Moreover, it changes every inequality in A to the * -C * form '<= '. * -C * * -C * Assumption: within a column in JA, the rows associated * -C * with nonzero entries are sorted in crescent order.* -C * * -C * called by: Subroutine PREPRO subroutines called: NONE* -C *********************************************************** -C -C ********************************************************** -C * Subroutine ONEONE, is called whenever a row has only * -C * one nonzero coefficient, moreover the corres- * -C * ponding column also has only one nonzero * -C * element. It results that either the problem * -C * is unbounded or infeasible or the value of * -C * the variable can be fixed. * -C * * -C * called by: subroutine PREPRO * -C * subroutines called: CHKFRE, UPTVFR * -C ********************************************************** -C -C ********************************************************** -C * Subroutine ONENZR, verifies for every nonzero element * -C * of the current row, which will be dropped, if * -C * it is the unique nonzero element in the corres- * -C * ponding column. If yes, the associated variable * -C * can be fixed or the problem is unbounded. * -C * * -C * called by: subr. PREPRO subroutines called: FIXVAR * -C * * -C * Remark: each row in JAT is sorted in crescent order of * -C * column index ( a result of the ARWISE subrout.)* -C ********************************************************** -C -C ********************************************************** -C * Subroutine DUALIZ, selects which problem to solve, * -C * * -C * * -C * Should work with just NCOLS, NROWS, also switches IBLA* -C * * -C * called by: Main subroutines called:None * -C ********************************************************** -C -C ********************************************************** -C * Subroutine UPTVFR, updates the array VRFREE if one * -C * finds a constraint which bounds a free * -C * variable. The variable will be no longer * -C * unrestricted and is deleted from VRFREE * -C * * -C * called by: subroutine PREPRO subroutines called:NONE * -C ********************************************************** -C -C ********************************************************** -C * Subroutine UPTCLD, updates arrays AT, JAT, IAT when * -C * a column is dropped. * -C * * -C * called by: Subroutine PREPRO subroutines called:NONE.* -C ********************************************************** -C -C ********************************************************** -C * SUBROUTINE UPTRWD, updates the arrays AT, JAT and IAT * -C * when a row is dropped. * -C * * -C * called by: subroutine PREPRO subroutines called: NONE* -C ********************************************************** -C -C ********************************************************** -C * Subroutine CHKRWS, checks if there are still variables* -C * in a constraint after a variable is fixed. If * -C * not, checks if the row is consistant and drop it* -C * if yes, otherwise the problem is infeasible. * -C * * -C * called by PREPRO subroutines called:DROPRW * -C ********************************************************** -C -C ********************************************************** -C * Subroutine CHKSGN, the parameter FRE returns the value* -C * zero if the given constraint does NOT have * -C * any free variable and all the nonzero * -C * coefficients have the same sign. Otherwise * -C * the returned value of FRE is POSITIVE. * -C * * -C * called by: subrout. PREPRO subroutines called:CHKFRE* -C ********************************************************** -C -C ********************************************************* -C * Subroutine CHKFRE, verifies if a given variable is or* -C * not unrestricted in sign. If yes, the * -C * parameter FRE returns the position occuppied * -C * by this variable in the array VRFREE. * -C * * -C * CALLED BY: sub. PREPRO,CHKSGN subrout. called: None * -C ********************************************************* -C -C ************************************************************************ -C * * -C * SWITCH * -C * * -C * SWITCH takes selected elements from the original * -C *data set and transfers them into the data * -C *structures required for the program. Once this -C *occurs, we process and then add pointers to * -C *each block. SWITCH also does portions of work previously found in AR- * -C * WISE * -C * * -C * Called by: main Calls: none * -C * Params: In, block number * -C * * -C ************************************************************************ -C * * -C * Pointers are initialized/handled as follows: * -C * rows: * -C * NROWS = total number of rows thus far, including preprocessing * -C * ITOTR = total number of rows read in from T data set * -C * * -C * Each is incremented in Switch, but only NROWS is inc./dec. when * -C * preprocessing occurs. Each "block" which is preprocessed goes * -C * from IBLA(.) to NROWS, and has NROWS-IBLA(.)+1 elements. * -C * * -C * Note that the order within columns is maintained, and each block * -C * starts with row 1, column 1. * -C * * -C ************************************************************************ -C -C Columns are "taken" from TIA, etc. in the order determined in SORTCL. -C This makes the non zero structure look like: -C -C AAAAA This way, we won't have to store zeroes in between the cols -C AAAAA of T. tsort(i) is the number of the ith column that has non -C TTT non zeroes in both A and T -C TTT -C -C !!!!!!!!!!!!!! NOTE NOTE NOTE !!!!!!!!!!!!!!!! -C -C This version (3/1/90) DOES NOT DUALIZE, unlike the previous version. -C -C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -C - SUBROUTINE SWITCH(Iblk) - INCLUDE (sizes.h) -C - CHARACTER*3 tRTYPE, RTYPE, BLANK1 - CHARACTER*8 tNAMER, tNAMEC, NAMERW, NAMECL, BLANK2 -C - DOUBLE PRECISION tA, tB, tC, VALFIX, VALLBD, ZEROD, - 1 A,B,C,ZCONST,X,Y,AT,taux -C - LOGICAL TEST,TEST2,TEST3 - INTEGER ISCEN,NAROWS,NACOLS,NWROWS,NWCOLS,IBLA,IBLAT,IBLAAT, - 1 IDIAG,IBLSZ,ITOTR,ITOTC,NJA,NJAT,IACLOF - INTEGER I, tNROWS, tNCOLS, IA, JA, ZEROI, NRWDRP, NFIX, NLBD, - 1 ORGIDX, VARFIX, VARLBD, NSTVAR, VRFREE, NFREE, - 1 VRCSGN, NCHSGN,tia, tja,iblk, nia, niat, h, I1,I2, - 1 tsort, nsort, NCOLS, NROWS, IAT, JAT, toridx -C - COMMON /debug/ idbg - COMMON /stoch1/ iscen, narows,nacols,nwrows,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /tCONST/ tNAMER(SIZEVEC), tRTYPE(SIZEVEC), tNAMEC(SIZEVEC) - COMMON /BLANK/ ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /tINPT1/ tA(SIZEMAT),tNROWS,tNCOLS,tIA(SIZEVEC), - 1 tJA(SIZEMAT),tB(SIZEVEC), tC(SIZEVEC), - 1 taux(SIZEVEC), tsort(SIZEBLK),NSORT, - 2 toridx(SIZEVEC),imxsrt - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) -C - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /INPT2/ B(SIZEVEC), C(SIZEMAT), ZCONST - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /FIXED/ VALFIX(SIZEBLK),VALLBD(SIZEBLK),VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK),VRFREE(SIZEVEC),VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C -C Test for A block it should be the first -C - IF(IBLK.EQ.1) THEN -C -C MOVE ROWS,COLUMNS -C - IBLA(1)=1 - NIA=NACOLS - NIAT=ZEROI - NJAT=ZEROI - NCOLS=NACOLS - NROWS=NAROWS - ITOTR=ZEROI -C -C - ITOTC=ZEROI - NJA=ZEROI - IA(1)=TIA(1) - IAT(1)=1 -C -C First, do columns not in TSORT, then do columns in TSORT -C - TEST=.TRUE. -C - 8 DO 22 I=1,NACOLS - TEST2=.FALSE. - DO 9 H=1,NSORT - IF(I.EQ.TSORT(H)) TEST2=.TRUE. - 9 CONTINUE - IF(TEST.AND.TEST2) GOTO 22 - IF(.NOT.TEST.AND..NOT.TEST2) GOTO 22 -C - ITOTC=ITOTC+1 - NAMECL(ITOTC)=tNAMEC(I) - ORGIDX(ITOTC)=I - C(ITOTC)=TC(I) - TEST3=.TRUE. -C - DO 10 J=TIA(I),TIA(I+1)-1 -C ....Only want elements with rows in A.... - IF(TJA(J).GT.NAROWS) GOTO 10 - TEST3=.FALSE. - NJA=NJA+1 - JA(NJA)=TJA(J) - A(NJA)=TA(J) - 10 CONTINUE - IF(TEST3) THEN - WRITE(6,*) 'Col with elements in T/A: ', - 1 'No. ',I, ' Name ',NAMECL(ITOTC) -C NJA=NJA+1 -C JA(NJA)=TJA(TIA(I)) -C WAS JA(NJA)=TJA(TIA(I)) -C A(NJA)=ZEROD - ENDIF - IA(ITOTC+1)=NJA+1 - 20 CONTINUE - 22 CONTINUE - IF(TEST) THEN - TEST=.FALSE. - GOTO 8 - ENDIF - IBLA(IBLK+1)=NIA+1 -C -C DO ROWS, AT -C - IBLAT(1)=1 -C -C objective -C - NAMERW(1)=TNAMER(1) - DO 30 I=1,NAROWS - NIAT=NIAT+1 - ITOTR=ITOTR+1 - NAMERW(I+1)=TNAMER(I+1) - RTYPE(I)=TRTYPE(I) -C Sign of b and A already changed in READetc. - B(I)=TB(I) -C - DO 28 I1=1,NACOLS - IEND= IA(I1+1) - 1 - DO 26 I2= IA(I1),IEND - if(idbg.ge.3) write(7,953) i2,ja(i2),j - IF (JA(I2) .EQ. I) THEN - NJAT=NJAT+1 - JAT(NJAT)= I1 - AT(NJAT)= A(I2) - GOTO 28 - ELSE - IF (JA(I2) .GT. I) GOTO 28 - ENDIF - 26 CONTINUE - 28 CONTINUE - IAT(NIAT+1)= NJAT+1 - 30 CONTINUE - rtype(narows+1)=trtype(narows+1) - IBLAT(2)=NIAT+1 -C - GOTO 999 - ENDIF -C -CTTTTTT -C IF 2 <= IBLK <= ISCEN+1 THEN WE HAVE A T MATRIX -c IF ISCEN+2 <= IBLK <= 2*ISCEN+1 THEN WE HAVE A W MATRIX -C - IF(IBLK.LE.ISCEN+1) THEN -C - ISTART=NAROWS+(IBLK-2)*NWROWS+1 - IEND=NAROWS+(IBLK-1)*NWROWS -C - DO 60 H=1,NSORT - I=TSORT(H) - NIA=NIA+1 - TEST=.TRUE. - DO 50 J=TIA(I),TIA(I+1)-1 -C ....Only want elements with rows iN T(IBLK) - - IF((TJA(J).LT.ISTART).OR.(TJA(J).GT.IEND)) GOTO 50 - TEST=.FALSE. - NJA=NJA+1 - JA(NJA)=TJA(J)-ISTART+1 - IF(JA(NJA).LE.0) THEN - WRITE(7,41) NJA,JA(NJA) - 41 FORMAT(' XXX - ERROR, JA(',I5,') = ',I5) - STOP - ENDIF - A(NJA)=TA(J) - 50 CONTINUE - IF(TEST) THEN - NJA=NJA+1 - JA(NJA)=1 - A(NJA)=ZEROD - ENDIF - IA(NIA+1)=NJA+1 - 60 CONTINUE - IBLA(IBLK+1)=NIA+1 -C -C DO AT, ROWS -C - DO 69 I=1,NWROWS -C - ITOTR=ITOTR+1 - NIAT=NIAT+1 - NROWS=NROWS+1 -C -C RTYPE is offset by one becoause of the objective row -C - NAMERW(NROWS+1)=TNAMER(ITOTR+1) - RTYPE(NROWS+1)=TRTYPE(ITOTR+1) -C /* keep one ahead for ISLACK */ -c rtype(nrows+1)=trtype(itotr+2) - B(NROWS)=TB(ITOTR) -C - DO 68 I1=IBLA(IBLK),IBLA(IBLK+1)-1 - IEND= IA(I1+1) - 1 - DO 66 I2= IA(I1),IEND - if(idbg.ge.3) write(7,953) i2,ja(i2),j - IF (JA(I2) .EQ. I) THEN - NJAT=NJAT+1 - JAT(NJAT)= I1-IBLA(IBLK)+1 - AT(NJAT)= A(I2) - GOTO 68 - ELSE - IF (JA(I2) .GT. I) GOTO 68 - ENDIF - 66 CONTINUE - 68 CONTINUE - IAT(NIAT+1)= NJAT+1 - 69 CONTINUE - IBLAT(IBLK+1)=NIAT+1 - GOTO 999 -C - ENDIF -C -CWWWWW So, it must be a W matrix..... -C -> No new rows need to be processed, only new columns -C - IF(IBLK.GT.2*ISCEN+1) GOTO 100 -C - ISTART=NAROWS+(IBLK-ISCEN-2)*NWROWS+1 - II1=NACOLS+(IBLK-ISCEN-2)*NWCOLS+1 - II2=NACOLS+(IBLK-ISCEN-1)*NWCOLS - DO 80 I=II1,II2 - ITOTC=ITOTC+1 - NCOLS=NCOLS+1 - NIA=NIA+1 - NAMECL(ncols)=tNAMEC(ITOTC) - ORGIDX(ncols)=i - C(ncols)=TC(ITOTC) -C -C WRITE(7,71) I,NIAt,NCOLS,ISTART -C 71 FORMAT(/,'I,NIAt,NCOLS,ISTART:',4(I5,2X)/) -C - DO 75 J=TIA(I),TIA(I+1)-1 - NJA=NJA+1 - JA(NJA)=TJA(J)-ISTART+1 - A(NJA)=TA(J) - 75 CONTINUE - IA(NIA+1)=NJA+1 - 80 CONTINUE - IBLA(IBLK+1)=NIA+1 -C -C AT(W) -C - DO 89 I=1,NWROWS -C - NIAT=NIAT+1 - DO 88 I1=IBLA(IBLK),IBLA(IBLK+1)-1 - IEND= IA(I1+1) - 1 - DO 86 I2= IA(I1),IEND - if(idbg.ge.3) write(7,953) i2,ja(i2),j - IF (JA(I2) .EQ. I) THEN - NJAT=NJAT+1 - JAT(NJAT)= I1-IBLA(IBLK)+1 - AT(NJAT)= A(I2) - GOTO 88 - ELSE - IF (JA(I2) .GT. I) GOTO 88 - ENDIF - 86 CONTINUE - 88 CONTINUE - IAT(NIAT+1)= NJAT+1 - 89 CONTINUE - IBLAT(IBLK+1)=NIAT+1 -C - GOTO 999 - 100 WRITE(6,110) IBLK - 110 FORMAT(' XXX - Error: Subroutine SWITCH: IBLK = ',I5) - STOP -C - 999 IF(IDBG.GE.1) THEN - WRITE(7,900) IBLK,ITOTR,ITOTC,ITOTR,NCOLS,NROWS,NJAt,NJA,NIA - 900 FORMAT(' SWITCH: IBLK,ITOTR,ITOTC,ITOTR,NCOLS,NROWS - 1, NJAt,NJA,NIA',9i5) - ENDIF - IF(IDBG.GE.3) THEN - DO 910 I=1,NCOLS - WRITE(7,950) I, IA(I),NAMECL(I) - DO 908 J=IA(I),IA(I+1)-1 - WRITE(7,951) J,JA(J),A(J) - 908 CONTINUE - 910 CONTINUE - DO 915 I=1,NROWS - WRITE(7,952) I, NAMERW(I),RTYPE(I) - 915 CONTINUE - 950 FORMAT (' COL,IAt,NAME:',2I5,A10) - 951 FORMAT (' J,JAt(J),A(J): ',8(I5,2X,I5,F5.2)) - 952 FORMAT (' I,NAMERW,RTYPE', 3(I5,2X,A10,2X,A3)) - 953 FORMAT (' ARWISE: I2,JAt(I2),J:',3(2X,I5)) -C - ENDIF -C - RETURN - END ##### sed > BQ_Decomp/wrapup.f <<'#####' 's/^-//' -C -C *************************************************************** -C * Subroutine GETSTX obtains the primal solution of the * -C * original problem. Essentially, it INCLUDEs to * -C * the solution obtained the varibles whose values * -C * were fixed as well as it updates the values of * -C * those variables whose bounds were changed during* -C * the solution procedure. * -C * * -C * called by: Main Program subroutines called: None * -C *************************************************************** -C - SUBROUTINE GETSTX(X,NCOLS,ZEROI,NROWS,Y,PROBLM) - INCLUDE (sizes.h) -C - CHARACTER*4 PROBLM -C - DOUBLE PRECISION VALFIX,VALLBD,AUX(SIZEVEC),X(SIZEVEC),Y(SIZEMAT) -C - INTEGER VARFIX, VARLBD, NFIX, NLBD, NSTVAR, ORGIDX, IEND, IBEG, - 1 NCOLS, ZEROI, VRFREE, NFREE, NROWS, NRWDRP, K1, K2, I, - 1 VRCSGN, NCHSGN -C - COMMON /FIXED/ VALFIX(SIZEBLK), VALLBD(SIZEBLK), VARFIX(SIZEBLK), - 1 VARLBD(SIZEBLK), VRFREE(SIZEVEC), VRCSGN(SIZEBLK), - 1 NFREE, NFIX, NLBD, NSTVAR, NRWDRP, NCHSGN, - 1 ORGIDX(SIZEVEC) -C - IF (PROBLM .EQ. 'DUAL') GOTO 40 - DO 10 I=1,NCOLS - AUX(I)=X(I) - 10 CONTINUE - DO 20 I=1,NROWS - X(I)=Y(I) - 20 CONTINUE - DO 30 I=1,NCOLS - Y(I)= AUX(I) - 30 CONTINUE - AUX(1)= NROWS - NROWS= NCOLS - NCOLS= AUX(1) - GOTO 90 -C - 40 DO 80 I=1,NFREE - IEND=NCOLS - K1= ZEROI - K2= ZEROI - DO 50 J=1,IEND - IF (VRFREE(I) .NE. ORGIDX(J)) GOTO 50 - IF (K1 .NE. ZEROI) GOTO 70 - K1= J - 50 CONTINUE - WRITE(6,60) VRFREE(I) - 60 FORMAT('The 2 non-negative components of',I4,' were not found.') - GOTO 180 -C - 70 K2= J - X(K1)= X(K1) + X(K2) -C................................. To compile 3/15/90 -C....... CALL DROPCL(K2) -C................................. - 80 CONTINUE -C - 90 IF (NFIX .EQ. ZEROI) GOTO 140 - IEND= NSTVAR - NFIX + 1 - IBEG= NCOLS -C -------------------------------------- -C REMARK: these components of X from IEND to IBEG are SLACK var. -C ---------------------------------------------------------- - DO 100 I= IBEG,IEND,-1 - AUX(I+NFIX)= X(I) - 100 CONTINUE -C - IEND= NSTVAR - NFIX - DO 110 I=1,IEND - AUX(ORGIDX(I))= X(I) - 110 CONTINUE -C - DO 120 I=1,NFIX -C write(6,444) i,varfix(i) -C 444 format('varfix(',i3,')= ',i3) - AUX(VARFIX(I))= VALFIX(I) - 120 CONTINUE -C - NCOLS= NCOLS+NFIX - DO 130 I=1,NCOLS - X(I)= AUX(I) - 130 CONTINUE -C - 140 IF (NLBD .EQ. ZEROI) GOTO 160 - DO 150 I=1,NLBD - X(VARLBD(I))= X(VARLBD(I)) + VALLBD(I) - 150 CONTINUE -C - 160 IF (NCHSGN .EQ. ZEROI) GOTO 180 - DO 170 I=1,NCHSGN - X(VRCSGN(I))= -1*X(VRCSGN(I)) - 170 CONTINUE -C - 180 RETURN - END -C -C ************************************************************* -C * * -C * SUBROUTINE PRINT : used to print out the final results, * -C * * -C ************************************************************* -C - SUBROUTINE PRINT(ZUPPER,ZLOWER,FEASBL,PROBLM,NFIX,NRWDRP, - 1 NSTFRE,PHASE,ARTFCL,BIGM) - INCLUDE (sizes.h) -C - CHARACTER*3 RTYPE - CHARACTER*4 PROBLM - CHARACTER*8 NAMERW, NAMECL, FEASBL -C - DOUBLE PRECISION A,B,C,Y,SLACK,ARTFCL,BIGM - DOUBLE PRECISION APRINT(20,20) - DOUBLE PRECISION DY,AAT,DIAG,LSTPRD,X,ZUPPER,ZLOWER,ZCONST -C - REAL PREPTM, READTM, STUPTM, PHS1TM, PHS2TM, ENDTIM -C - INTEGER NROWS,NCOLS,IA,JA,IAAT,JAAT,PTRPRD,LSTDST - INTEGER I,J,IEND,LAST, NRWDRP, NFIX, PHASE, NSTFRE - INTEGER ITERA1, ITERA2 -C - common /stoch1/ iscen, narows,nacols,nwrols,nwcols, ibla(MAXBLK), - 1 iblat(MAXBLK),iblaat(MAXBLK),idiag(MAXBLK), - 1 iblsz(MAXBLK,2), itotr, itotc, nja,nia,niat,njat, - 1 iaclof - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /INPT2/ B(SIZEVEC), C(SIZEMAT), ZCONST - COMMON /INPT3/ Y(SIZEVEC), X(SIZEVEC) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - COMMON /OUTPRO/ LSTPRD(SIZEVEC),PTRPRD(SIZEMAT),LSTDST(SIZEMAT) - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, STUPTM, - 1 PHS1TM, PHS2TM, MXITNS - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) -C - IF (FEASBL .EQ. 'FEASIBLE') GOTO 210 - IF ((FEASBL .EQ. 'UNBOUNDD') .AND. (PHASE .EQ. 2)) THEN - WRITE(6,780) - GOTO 210 - ENDIF - IF ((FEASBL .EQ. 'UNBOUNDD') .AND. (PHASE .EQ. 1)) THEN - WRITE(6,790) - GOTO 210 - ENDIF - IF (FEASBL .EQ. 'STUCK-1 ') THEN - WRITE(6,800) - GOTO 210 - ENDIF - IF (FEASBL .EQ. 'STUCK-2 ') THEN - WRITE(6,850) ARTFCL - GOTO 210 - ENDIF - WRITE(6,770) FEASBL - GOTO 1000 -C ------------------------------------------------------------ -C IF (PROBLM .EQ. 'PRML') THEN -C IEND1= NROWS -C IEND2= NCOLS -C ELSE -C IEND1=NCOLS -C IEND2=NROWS -C ENDIF -C WRITE(*,510) NAMERW(1),(C(I),I=1,IEND1) -C -C DO 130 I1=1,IEND2 -C DO 120 I2=1,IEND1 -C APRINT(I1,I2)=0.0D0 -C 120 CONTINUE -C 130 CONTINUE -C DO 170 I1=1,IEND1 -C DO 160 I2=IA(I1),IA(I1+1)-1 -C LAST=1 -C 140 IF (LAST .EQ. JA(I2)) GOTO 150 -C LAST= LAST+1 -C GOTO 140 -C 150 APRINT(LAST,I1)=A(I2) -C 160 CONTINUE -C 170 CONTINUE -C -C DO 180 I=1,IEND2 -C WRITE(*,520) NAMERW(I+1),(APRINT(I,J),J=1,IEND1),B(I) -C 180 CONTINUE -C -C WRITE(6,601) (SLACK(I),I=1,IEND1) -C WRITE(6,760) -C WRITE(6,690) -C ----------------------------------------------------------- -C 10 WRITE(6,690) -C DO 190 I=1,NROWS,3 -C IF ((I+2) .GT. NROWS) GOTO 185 -C WRITE(6,710) I,Y(I),I+1,Y(I+1),I+2,Y(I+2) -C GOTO 190 -C 185 WRITE(6,710) (J,Y(J),J=I,NROWS) -C 190 CONTINUE -C WRITE(6,760) -C WRITE(6,700) -C DO 210 I=1,NCOLS,3 -C IF ((I+2) .GT. NCOLS) GOTO 200 -C WRITE(6,720) I,X(I),I+1,X(I+1),I+2,X(I+2) -C GOTO 210 -C 200 WRITE(6,720) (J,X(J),J=I,NCOLS) - 210 CONTINUE - WRITE(6,760) - WRITE(6,730) ZLOWER - WRITE(6,740) ZUPPER - IF (ARTFCL .LE. 0.0D0) GOTO 220 - WRITE(6,640) - IF (PROBLM .EQ. 'DUAL') THEN - WRITE(6,730) ZLOWER+ARTFCL*BIGM - WRITE(6,740) ZUPPER - ELSE - WRITE(6,730) ZLOWER - WRITE(6,740) ZUPPER-ARTFCL*BIGM - ENDIF - WRITE(6,820) BIGM - WRITE(6,830) ARTFCL - WRITE(6,760) - 220 WRITE(6,480) NROWS+NRWDRP - WRITE(6,490) NCOLS - WRITE(6,500) ZCONST - WRITE(6,530) NRWDRP - WRITE(6,540) NFIX - WRITE(6,570) NSTFRE - WRITE(6,550) PROBLM - WRITE(6,560) PHASE -C WRITE(6,580) PREPTM -C WRITE(6,590) READTM -C WRITE(6,600) STUPTM - WRITE(6,610) -C WRITE(6,620) PHS1TM, PHS2TM - WRITE(6,750) ITERA1, ITERA2 - ENDTIM= DTIME(0) -C WRITE(6,650) ENDTIM -C WRITE(6,630) READTM+PREPTM+STUPTM+PHS1TM+PHS2TM+ENDTIM -C - 480 FORMAT(10X,'# of ORIGINAL constraints ..',I14) - 490 FORMAT(10X,'# of ORIGINAL columns ......',I14) - 500 FORMAT(10X,'constant in Z ..............',F14.1) - 510 FORMAT(A8,13F5.1//) - 520 FORMAT(A8,14F5.1) - 530 FORMAT(10X,'# of constraints dropped ...',I14) - 540 FORMAT(10X,'# of variables fixed .......',I14) - 550 FORMAT(10X,'Problem solved .............',A14) - 560 FORMAT(10X,'Stopped in ................. PHASE',I2) - 570 FORMAT(10X,'# of free variables ........',I14) - 580 FORMAT(10X,'Preprocessing time .........',F14.5) - 590 FORMAT(10X,'Reading input data .........',F14.5) - 600 FORMAT(10X,'Setting up time ............',F14.5,/) - 601 FORMAT('DUAL SLK',13F5.1) - 610 FORMAT(30X,'Phase 1',8X,'Phase 2') - 620 FORMAT(10X,'Time .......',5X,F10.5,9X,F10.5) - 630 FORMAT(10X,'Total time .................',F14.5,/) - 640 FORMAT(10X,'Bounds without the ARTIFICIAL variable, ',/) - 650 FORMAT(10X,'Final Calc. & Printing .....',F14.5) - 690 FORMAT(//,3X,'DUAL SOLUTION') - 700 FORMAT(/,3X,'PRIMAL SOLUTION') - 710 FORMAT(/,3(3X,'Y(',I3,')= ',E11.4)) - 720 FORMAT(/,3(3X,'X(',I3,')= ',E11.4)) - 730 FORMAT(10X,'Lower Bound on Z ...... ZLOWER= ',E15.9) - 740 FORMAT(10X,'Upper Bound on Z ...... ZUPPER= ',E15.9,/) - 750 FORMAT(10X,'Iterations..',10X,I5,10X,I5) - 760 FORMAT(' ') - 770 FORMAT('The constraint (column) ',A8, - 1' makes this problem INFEASIBLE (UNBOUNDED).') - 780 FORMAT('The problem solved is UNBOUNDED. The current solution is') - 790 FORMAT('WARNING. Can this problem be UNBOUNDED and INFEASIBLE?') - 800 FORMAT(/,10X,'WARNING. The step size has reached its upper', - 1' bound.',/,10X,' Either the algorithm is trapped near the', - 1' boundary',/,10X,'or the problem being solved is unbounded.', - 1/,10X,'The current solution is:'//) - 820 FORMAT(/,10X,'BIGM = ',E15.9) - 830 FORMAT(10X,'ARTIFICIAL = ',E15.9,/) - 850 FORMAT(/,10X,'WARNING. The value of the ARTIFICIAL variable is ', - 1/,10X,E14.7,' and decreases very slowly. Either',/,10X, - 1'the problem being solved is INFEASIBLE or the',/,10X, - 1'the algorithm is trapped near the boundary.Maybe',/,10X, - 1'you should try a different starting point or to',/,10X, - 1'solve the dual problem. The current solution is:',//) -C - 1000 RETURN - END #####