Shell the following script in an empty directory to create the for-files associated with the paper: "Efficient solution of two stage stochastic linear programs using interior point methods" by J.R. Birge and D.F. Holmes. mkdir DualAffCode mkdir BQ_Decomp sed > README <<'#####' 's/^-//' -The following two files contain the codes used to generate -results for the paper -"Efficient Solution of Two Stage Stochastic Linear Programs using -Interior Point Methods", by J.R. Birge, and D.F. Holmes: - -File Notes -DualAffCode directory that contains source codes implementing the - dual affine scaling scheme (with Schur Complement) - used in the paper. Usage requires two packages: - SPARSPAK and LINPACK (see references in paper). - -BQ_Decomp directory that contains source codes for an implementation - of the Birge-Qi decomposition used in the paper. Requires - SPARSPAK and LINPACK. - - Credits: Both codes use routines originally developed by Jose Arantes - (Univ. of Cincinnati). - -Derek Holmes -Industrial and Ops. Engineering -Univ. of Michigan -Questions? email dholmes@engin.umich.edu - - ##### sed > DualAffCode/bldupd.f <<'#####' 's/^-//' -C *********************************************************** -C * Subroutine COMPCT: sorts IAUX in crescent order and * -C * eliminates duplicated entries. * -C *********************************************************** -C - SUBROUTINE COMPCT(IAUX,TOTAL) - INCLUDE (common.h) -C - INTEGER IAUX(MXRWS), TOTAL, TTSEEN -C - TTSEEN= 1 - 10 IF (TTSEEN .GE. TOTAL) GOTO 40 - MIN= IAUX(TTSEEN) - IDXMIN= TTSEEN -C - DO 30 I= TTSEEN+1, TOTAL - 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 * Subroutine INIMTR: Initializes the values of the matrix * -C * ADk2At which is stored through the * -C * arrays DIAG, AAT, IAAT, JAAT . * -C * * -C * called by: Main Program subroutines called: None * -C *************************************************************** -C - SUBROUTINE INIMTR - INCLUDE (common.h) -C - INTEGER flag -C - CALL FINDDG - K=1 - flag=0 - AAT(K)=ZEROD - 10 DO 90 I1= 1, NROWS-1 - IAAT(I1)=K - DO 80 I2= (I1+1),NROWS - DO 70 L= 1,NSPCOLS - DO 20 I3= IA(L),IA(L+1)-2 - IF (JA(I3) .EQ. I1) GOTO 30 - 20 CONTINUE - GOTO 70 - 30 M1=I3 - DO 40 I3= (M1+1),IA(L+1)-1 - IF (JA(I3) .EQ. I2) GOTO 50 - 40 CONTINUE - GOTO 70 - 50 M2=I3 - LL=MPSP2C(L) - AAT(K)=AAT(K)+A(M1)*A(M2)*D(LL) - flag=1 - 70 CONTINUE - IF (flag .eq. 0) GOTO 80 - flag=0 - JAAT(K)= I2 - K= K+1 - AAT(K)= ZEROD - 80 CONTINUE - 90 CONTINUE - IAAT(NROWS)= K -C - DO 100 I1=1,NROWS - DIAG(I1)= ZEROD - 100 CONTINUE - DO 120 I1=1,NSPCOLS - I2=MPSP2C(I1) - DO 110 J2=IA(I1),IA(I1+1)-1 - DIAG(JA(J2))=DIAG(JA(J2))+A(J2)*A(J2)*D(I2) - 110 CONTINUE - 120 CONTINUE - nzerop=zeroi - do 130 i1=1,nrows - if(diag(i1).le.epsilo) then - diag(i1)=EPSILO - nzerop=nzerop+1 - endif - 130 continue -C - RETURN - END -C -C ************************************************************ -C * * -C * Subroutine OUTPRD: One computes every nonzero * -C * OUTER PRODUCT A(i,l) x A(j,l), * -C * once in the beginning of the LP procedure. With * -C * these values, this subroutine fills in the arrays * -C * PTRPRD, LSTPRD and LSTDST, which will help at each * -C * iteration, to buils the matrix Bk= ADk2At. * -C * * -C * called by: MAIN program * -C * * -C * subroutines called: None * -C ************************************************************ -C - SUBROUTINE OUTPRD(jerror) - INCLUDE (common.h) -C -------------------------------------------------------------- -C PTRPRD: contains pointers to the first nonzero element of -C the outer-product list involving elements in each -C column of A -C LSTPRD: contains the computed outer-products sorted in -C increasing order of the corresponding column index l -C LSTDST: for each entry in LSTPRD, LSTDST contains a pointer -C to the location of the element Bk(i,j) in array AAT -C ---------------------------------------------------------------- - L=1 - DO 130 J=1,NSPCOLS - PTRPRD(J)= L - IAUX2= IA(J+1) - 2 - IAUX1= IA(J+1) - 1 - IF (IAUX2 .LT. IA(J)) GOTO 130 - DO 120 I= IA(J),IAUX2 - DO 110 K=(I+1),IAUX1 - LSTPRD(L)= A(I)*A(K) - DO 10 M=IAAT(JA(I)),(IAAT(JA(I)+1)-1) - IF (JAAT(M) .EQ. JA(K)) GOTO 30 - 10 CONTINUE - WRITE(6,920) - write(6,915) K,JA(K) - WRITE(6,*) ' Trying to find JA(',K,') = ',JA(K) - write(6,*) ' COLumn: ',j,' Name: ',namecl(j) - write(6,*) ' jaat(iaat(ja(i))))-jjat(iaat(.+1)-1)' - write(6,*) (n,jaat(n),' : ',n=iaat(ja(i)), - 1 iaat(ja(i)+1)-1 ) - error=1 - GOTO 140 - 30 LSTDST(L)= M - L= L+1 - 110 CONTINUE - 120 CONTINUE - 130 CONTINUE - PTRPRD(NSPCOLS+1)=L -C - 140 RETURN - 915 format('JA(',I3,')= ',I4) - 920 FORMAT('ERROR: NO COLUMN IN Bk ASSOCIATES TO OUTER PRD') - END -C - -C -C *************************************************************** -C * * -C * Subroutine UPDMTR : This subroutine, at each iteration * -C * updates the matrix Bk= ADkDkAt * -C * However the only elements of Dk changing are those * -C * whose relative change in this iteration is greater * -C * then EPSUPD. * -C * * -C * called by: MAIN PROGRAM Subroutines called: None * -C *************************************************************** -C - SUBROUTINE UPDMTR - INCLUDE (common.h) -C ---------------------------------------------- Only update sparse - DO 30 I=1,NSPCOLS - MAP = MPSP2C(I) - AUX1= DBONE/SLACK( MAP ) - AUX2= DABS(AUX1*AUX1-D(MAP))/DABS(D(MAP)) - IF (AUX2 .LT. EPSUPD) GOTO 30 - AUX3= (AUX1*AUX1) - (D(MAP)) - DO 10 J=IA(I), IA(I+1)-1 - DIAG(JA(J))=DIAG(JA(J))+A(J)*A(J)*AUX3 - 10 CONTINUE - DO 20 J= PTRPRD(I), PTRPRD(I+1)-1 - AAT(LSTDST(J))=AAT(LSTDST(J))+LSTPRD(J)*AUX3 - 20 CONTINUE - 30 CONTINUE - CALL FINDDG -C - RETURN - END -C -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(COEFSK,IROW) - INCLUDE (common.h) -C - NSPCOLS= NSPCOLS+1 - IA(NSPCOLS+1)= IA(NSPCOLS)+1 - JA(IA(NSPCOLS))= IROW - A(IA(NSPCOLS))= COEFSK - C(NSPCOLS)= ZEROD - COLSF(NSPCOLS)=DBONE -C - DO 10 J= IROW+1, NROWS+1 - IAT(J)= IAT(J)+1 - 10 CONTINUE -C - DO 20 J= IAT(NROWS+1)-2, IAT(IROW+1)-1, -1 - JAT(J+1)= JAT(J) - AT(J+1)= AT(J) - 20 CONTINUE -C - JAT(IAT(IROW+1)-1)= NSPCOLS - AT(IAT(IROW+1)-1)= COEFSK -C - RETURN - END -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 - SUBROUTINE UPTCLD(ICLDROP) - INCLUDE (common.h) -C - NXTJAT= ZEROI - DO 30 I=1,NROWS - IAT(I)= NXTJAT+1 - DO 20 J= IAT(I), IAT(I+1)-1 - IF (JAT(J) .EQ. ICLDROP) GOTO 20 - IF (JAT(J) .LT. ICLDROP) GOTO 10 - JAT(J)= JAT(J)-1 - 10 NXTJAT= NXTJAT+1 - JAT(NXTJAT)= JAT(J) - AT(NXTJAT)= AT(J) - 20 CONTINUE - 30 CONTINUE - IAT(NROWS+1)= NXTJAT+1 -C - RETURN - END -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 - SUBROUTINE UPTRWD(IRWDROP) - INCLUDE (common.h) -C - INTEGER TOTNZC -C - TOTNZC= IAT(IRWDROP+1) - IAT(IRWDROP) - IF (IRWDROP .EQ. (NROWS+1)) GOTO 30 - DO 10 I=IRWDROP,NROWS - IAT(I+1)=IAT(I+2)-TOTNZC - 10 CONTINUE -C - DO 20 I= IAT(IRWDROP), IAT(NROWS+1)-1 - JAT(I)= JAT(I+TOTNZC) - AT(I)= AT(I+TOTNZC) - 20 CONTINUE -C - 30 DO 40 I= IAT(NROWS+1), IAT(NROWS+2) - JAT(I)= ZEROI - AT(I)= ZEROD - 40 CONTINUE - IAT(NROWS+2)= ZEROI -C - RETURN - END -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 - SUBROUTINE UPTVFR(IDX) - INCLUDE (common.h) -C - NFREE= NFREE-1 - IF (NFREE .LE. ZEROI) GOTO 20 - DO 10 I=IDX,NFREE - VRFREE(I)= VRFREE(I+1) - 10 CONTINUE -C - 20 RETURN - END ##### sed > DualAffCode/common.h <<'#####' 's/^-//' - PARAMETER ( MXNZS = 80000 ) - PARAMETER ( MXRWS = 6000 ) - PARAMETER ( MXCLS = 8000 ) - PARAMETER ( MXSNZ = 150000 ) - PARAMETER ( MXMSC = 500 ) - PARAMETER ( MXDNCO = 40 ) - PARAMETER ( MXDNNZ = 5000 ) - PARAMETER ( ISUBY = 1 ) - PARAMETER ( ISUBN = 0 ) - PARAMETER ( ISUBP = -1 ) - PARAMETER ( IDUALA = 1 ) - PARAMETER ( IPRDUA = 2 ) - PARAMETER ( MXFLDS = 2 ) - IMPLICIT DOUBLE PRECISION (A-H,P-Z) - IMPLICIT INTEGER (I-O) - CHARACTER*3 BLANK1,RTYPE - CHARACTER*4 PROBLM - CHARACTER*8 BLANK2, NAMERW, NAMECL, FEASBL - REAL PREPTM, READTM, STUPTM, PHS1TM, PHS2TM - INTEGER VRCSGN, VARFIX, VARLBD, PHASE, ZEROI, VRFREE, PTRPRD - DOUBLE PRECISION LSTPRD - LOGICAL KSPLIT, KSCALE, KSCHUR, IDFLAG -C - COMMON /GLOBL/ RCDNO(4), IDEB(10), ITIM(30), IALG, KSCALE, - 1 KSCHUR, KSPLIT - COMMON /SIZES/ NROWS, NSPCOLS, NDNCOLS, NCOLS, NZEROP, - 1 NFIX,NFREE,NLBD,NRWDRP,NSTVAR, NCHSGN - COMMON /BLANK/ DBONE, ZEROD,ZEROI,BLANK1,BLANK2 - COMMON /CONST/ EPSSTP, EPSPH1, EPSILO, EPSUPD, GAMMA, UBOUND, - 1 UBSTEP, SLKMIN, BIGM, BIGMIU, EPSSTY, MAXITS -C - COMMON /PDATA/ A(MXNZS), IA(MXCLS),JA(MXNZS), B(MXRWS), - 1 C(MXCLS), ZCONST - COMMON /CURPT/ Y(MXRWS), X(MXCLS), SLACK(MXCLS) - COMMON /ROWISE/ AT(MXNZS), IAT(MXCLS), JAT(MXNZS) - COMMON /ADK2AT/ AAT(MXNZS),DIAG(MXRWS),IAAT(MXRWS),JAAT(MXNZS) - COMMON /OUTPRO/ LSTPRD(MXNZS),PTRPRD(MXRWS),LSTDST(MXNZS) - COMMON /BODY/ DY(MXRWS), DV(MXCLS), D(MXCLS), DROOT(MXCLS),STEP -C - COMMON /SPAUSR/ MSGLVA,IERRA,MAXSA,NEQNS - COMMON /STATUS/ ZOLD, YOLD, ZUPPER, ZLOWER, ZACTLOW, ITERA1, - 1 ITERA2, PHASE, FEASBL - COMMON /CONSTR/ NAMERW(MXRWS), RTYPE(MXRWS), NAMECL(MXCLS), - 1 PROBLM - COMMON /FIXED/ VALFIX(MXMSC), VALLBD(MXMSC), VARFIX(MXMSC), - 1 VARLBD(MXMSC), VRFREE(MXMSC), VRCSGN(MXMSC), - 1 ORGIDX(MXCLS) -C - COMMON /SCALDT/ ROWSF(MXRWS), COLSF(MXCLS) - COMMON /SCHUR/ AD(MXDNNZ), ATD(MXDNNZ), IAD(MXDNCO), JAD(MXDNNZ), - 1 IATD(MXDNCO), JATD(MXDNNZ), IDFLAG(MXCLS), - 1 MPSP2C(MXCLS), MPDN2C(MXDNCO), IDNTHR ##### sed > DualAffCode/debug.f <<'#####' 's/^-//' -C DEBUG(1) : A, AT -C DEBUG(2) : Names -C DEBUG(3) : AAT -C DEBUG(4) : DIAG, SLACK -C DEBUG(5) : Y -C DEBUG(6) : STEPSZ, DY -C DEBUG(7) : DV -C DEBUG(8) : SCale - - SUBROUTINE PRTDEB(IOPT) - INCLUDE (common.h) -C - write(7,*) ' DEBUG ==================================== ',iopt -C - If(iopt.eq.1) then - write(7,*) ' NROWS,NCOLS, : NSPCOLS, NDNCOLS' - write(7,*) nrows, ncols, nspcols, ndncols - write(7,*) 'NFIX,NFREE,NLBD,NRWDRP,NSTVAR, NCHSG' - write(7,*) NFIX,NFREE,NLBD,NRWDRP,NSTVAR, NCHSG - write(7,*) '....................................... IA , JA, A' - do 770 ii=1,nspcols - write(7,8031) II, IA(II), MPSP2C(II) - write(7,801) (J,JA(J),J=IA(II),IA(II+1)-1) - write(7,802) (J,A(J),J=IA(II),IA(II+1)-1) - 770 continue - write(7,*) '....................................... IAD,JAD,AD' - do 771 ii=1,ndncols - write(7,8032) II, IAD(II), MPDN2C(II) - write(7,801) (J,JAD(J),J=IAD(II),IAD(II+1)-1) - write(7,802) (J,AD(J),J=IAD(II),IAD(II+1)-1) - 771 continue - write(7,*) '..................................... IAT, JAT, AT' - if(iat(1).ne.0) then - do 7711 ii=1,nrows - write(7,804) II, IAT(II) - write(7,801) (J,JAT(J),J=IAT(II),IAT(II+1)-1) - write(7,802) (J,AT(J),J=IAT(II),IAT(II+1)-1) - 7711 continue - write(7,*) ' ..................................... IATD, JATD' - do 7712 ii=1,nrows - write(7,804) II, IATD(II) - write(7,801) (J,JATD(J),J=IATD(II),IATD(II+1)-1) - write(7,802) (J,ATD(J),J=IATD(II),IATD(II+1)-1) - 7712 continue - - endif - endif -C - if(iopt.eq.2) then - do 772 i=1,nrows+1 - write(7,*) i, b(i), namerw(i), rtype(i) - 772 continue - do 773 i=1,nspcols+ndncols - write(7,*) i,c(i), namecl(i) - 773 continue - endif -C - if(iopt.eq.3) then - write(7,*) ' ............. IAAT, JAAT, AAT ........' - write(7,810) (I,IAAT(I),I=1,nrows) - do 774 i=1,nrows - write(7,811) (J,JAAT(J), J=IAAT(I),IAAT(I+1)-1 ) - 774 continue - do 776 i=1,nrows - write(7,812) (J,AAT(J), J=IAAT(I),IAAT(I+1)-1 ) - 776 continue - endif -C - if(iopt.eq.4) then - write(7,820) (I,DIAG(I),I=1,NSPCOLS+NDNCOLS) - write(7,821) (I,SLACK(I),I=1,NSPCOLS+NDNCOLS) - endif -C - if(iopt.eq.5) then - write(7,830) (I,Y(I),I=1,NROWS) - write(7,*) ' Y(NROWS+1) = ARTFCL = ',Y(NROWS+1) - endif -C - if(iopt.eq.6) then - write(7,*) ' STEPSIZE: ',STEP - write(7,831) (I,DY(I),I=1,NROWS+1) - endif -C - if(iopt.eq.7) then - write(7,832) (I,DV(I),I=1,NSPCOLS+NDNCOLS) - endif -C - if(iopt.eq.8) then - write(7,840) (I,ROWSF(I),I=1,NROWS) - write(7,841) (I,COLSF(I),I=1,NSPCOLS+NDNCOLS) - endif -C - 801 format(6(' (',I3,') ',I5)) - 802 format(6(' (',I3,') ',E15.9)) - 8031 format(' IA(',I5,'): ',I5, ' MPSP2C: ',I5) - 8032 format('IAD(',I5,'): ',I5, ' MPDN2C: ',I5) - 804 format(' IAT(',I5,'): ',I5) - 805 format(6(' b(',I3,') ',E15.9)) - 806 format(6(' c(',I3,') ',E15.9)) - 810 format(6(' IT(',I3,'): ',I4)) - 811 format(6(' JT(',I5,'): ',I5)) - 812 format(6(' AAT(',I3,') ',E15.9)) - 820 format(6(' DIA(',I3,') ',E15.9)) - 821 format(6(' SLK(',I3,') ',E15.9)) - 830 format(6(' Y(',I3,') ',E15.9)) - 831 format(6(' DY(',I3,') ',E15.9)) - 832 format(6(' DV(',I3,') ',E15.9)) - 840 format(6(' RSF(',I3,') ',E15.9)) - 841 format(6(' CSF(',I3,') ',E15.9)) -C - RETURN - END ##### sed > DualAffCode/finddy.f <<'#####' 's/^-//' -C *************************************************************************** -C * Subroutine FINDDY - Finds the direction of improvement using sparspack -C * and routines from linpack. -C * By dfh - 6/16/92 includes sparspak stuff, IMPDRT, and FINDV -C * To find phase I search direction: (1) Solve ADA^T w = ADe for w, -C * (2) Solve ADA^T u = b for u, (3) Let dy(a)=(-M+ADe.u)/(e.De - (ADe)^T w -C * Here, v=ADe. (4) Let dy=u+w*dy(a) -C **************************************************************************** - SUBROUTINE FINDDY(S) - INCLUDE (common.h) - DIMENSION S(MXSNZ), VV(MXRWS), WW(MXRWS) -C - CALL ADASLV(S,DY,B) - IF (PHASE .EQ. 1) THEN -C ............................................ FINDV: Calc. AD^2e - DO 72 I=1,NROWS - VV(I)= ZEROD - 72 CONTINUE - CALL FINDDG -C ....................................... Sparse Portion - DO 76 ICOL= 1,NSPCOLS - DO 74 J= IA(ICOL), IA(ICOL+1)-1 - VV(JA(J))= (A(J)*D( MPSP2C(ICOL)) )+VV(JA(J)) - 74 CONTINUE - 76 CONTINUE -C ....................................... Dense Portion - DO 86 ICOL= 1,NDNCOLS - DO 84 J= IAD(ICOL), IAD(ICOL+1)-1 - VV(JAD(J))= (AD(J)*D( MPDN2C(ICOL) ))+VV(JAD(J)) - 84 CONTINUE - 86 CONTINUE -C ................................................ IMPDRT -C ....................................... Solve ADA w = v - CALL ADASLV(S,WW,VV) - DVW = DDOT(NROWS,VV,1,WW,1) - DVDY = DDOT(NROWS,VV,1,DY,1) - DDD = DDOT(NCOLS,DROOT,1,DROOT,1) - DY(NROWS+1)= (B(NROWS+1)+DVDY)/(DDD - DVW) - CALL DAXPY(NROWS,DY(NROWS+1),WW,1,DY,1) -C - ENDIF - 901 format(2(I3,2x,E20.8)) - RETURN - END -C ***************************************************************************** -C * Subroutine ADASLV - Solves systems of the form ADA^T u = v -C * By dfh 6/20/92 -C * ........ Includes Schur Complement: If schur complement not used, -C * then all that's necessary is to use the first sparse colve. -C Routine: (1) Let V = Ad.Dd (2) Solve [I + Vt (LLt)^{-1} V ] delta= -C -Vt (LLt)^{-1} b. -C ******************************************************************************* - SUBROUTINE ADASLV(S,u,v) - INCLUDE (common.h) - DIMENSION S(MXSNZ), U(MXRWS), V(MXRWS), DLLTV(MXRWS,MXDNCO) - DIMENSION DSCH(MXDNCO,MXDNCO), DELTA(MXDNCO) - DIMENSION KVPT(MXDNCO), Z(MXDNCO) - LOGICAL DEBUGSPK - DEBUGSPK=.FALSE. -C ................................................ Input Numerical Values - IF(DEBUGSPK) MSGLVA=6 - DO 50 I=1, NROWS-1 - DO 40 J= IAAT(I),IAAT(I+1)-1 - CALL INAIJ5(JAAT(J), I, AAT(J), S) - 40 CONTINUE - 50 CONTINUE - DO 60 I=1,NROWS - CALL INAIJ5(I,I,DIAG(I),S) - 60 CONTINUE -C - DO 65 I=1,NROWS - CALL INBI(I,V(I),S) - 65 CONTINUE -C ................................................ Factor and Solve - IF(.NOT.DEBUGSPK) MSGLVA=1 - CALL SOLVE5(S) - CALL ESTCOND(RCDNO(1),S) - DO 70 I=1,NROWS - U(I)= S(I) - 70 CONTINUE -C - IF( KSCHUR ) THEN -C .......................................... Zero out matrix - DO 110 I1=1,NDNCOLS - DO 100 J1=1,NROWS - DLLTV(J1,I1)=ZEROD - 100 CONTINUE - 110 CONTINUE - DO 180 I1=1,NDNCOLS -C .................................. Solve LLt x=V(.,I1) & save - DO 120 J1=IAD(I1),IAD(I1+1)-1 - CALL INBI(JAD(J1), AD(J1)*DROOT( MPDN2C(I1) ), S) - 120 CONTINUE - CALL SOLVE5(S) - DO 130 J1=1, NROWS - DLLTV(J1,I1)=S(J1) - 130 CONTINUE -d do 131 k1=1,ndncols -d write(7,'(4(1x,E15.8))') (dlltv(k2,k1),k2=1,nrows) -d131 continue -C .................................. Calc. Vt (LLt)^{-1} V +I - DO 150 I2=1,NDNCOLS - DTEMP=ZEROD - DO 140 J1= IAD(I2), IAD(I2+1)-1 - DTEMP=DTEMP+AD(J1)*S( JAD(J1) ) - 140 CONTINUE - DSCH(I1,I2) = DTEMP * DROOT( MPDN2C(I2) ) - 150 CONTINUE - DSCH(I1,I1)=DSCH(I1,I1)+DBONE -d write(7,*) ' SCHUR complemenet........' -d do 151 k1=1,ndncols -d write(7,'(4(1x,E15.8))') (dsch(k2,k1),k2=1,ndncols) -d151 continue -C .................................. Calc. -Vt (LLt)^{-1} b - DELTA(I1)=ZEROD - DO 160 J1= IAD(I1), IAD(I1+1)-1 - DELTA(I1)=DELTA(I1) + AD(J1)*U( JAD(J1) ) - 160 CONTINUE - DELTA(I1)=-DBONE*DELTA(I1)* DROOT( MPDN2C(I1) ) - 180 CONTINUE -d write(7,*) ' RHS: ' -d write(7,'(5(I3,1X,E20.12,1X))') (i,delta(i),i=1,ndncols) -C ..................................... Solve DSCH delta = RHS - MAXTEMP=MXDNCO - CALL DGECO(DSCH, MAXTEMP, NDNCOLS, KPVT, RCDNO(2), Z) - CALL DGESL(DSCH, MAXTEMP, NDNCOLS, KPVT, DELTA, ZEROI) -d write(7,*) ' DELTA: REAL.....' -d write(7,'(5(I3,1X,E20.12,1X))') (i,delta(i),i=1,ndncols) -C ..................................... Form u=u+(LLt.V).delta - DO 200 I1=1,NDNCOLS - DO 190 J1=1,NROWS - U(J1) = U(J1) + DLLTV(J1,I1)*DELTA(I1) - 190 CONTINUE - 200 CONTINUE - ENDIF -C ................................................... D O N E - RETURN - END -C ***************************************************************************** -C * Subroutine INITSP - Sets up the relevant data into sparspak -C * By dfh 6/16/92 -C ******************************************************************************* - SUBROUTINE INITSP(S) - INCLUDE (common.h) - DIMENSION S(MXSNZ) - LOGICAL DEBUGSPK - DEBUGSPK=.FALSE. -C -------------------- -C Input Structure -C -------------------- - IF(DEBUGSPK) MSGLVA=6 - IF(.NOT.DEBUGSPK) MSGLVA=1 - write(6,*) ' No. nonzeroes in aat ',iaat(nrows)-1 - dens=dfloat(iaat(nrows)-1)/dfloat(nrows*ncols) - write(6,*) ' Density: ', dens - ITIM(5)=MCLOCK() - CALL IJBEGN - DO 20 I=1, NROWS-1 - DO 10 J= IAAT(I), IAAT(I+1)-1 - CALL INIJ(JAAT(J),I,S) - CALL INIJ(I,I,S) - 10 CONTINUE - 20 CONTINUE - CALL INIJ(NROWS,NROWS,S) - CALL IJEND(S) - IF(.NOT.DEBUGSPK) MSGLVA=1 -C ---------------------------------------- -C Determine Symmetric Ordering -C ( B5= Minimum Ordering Degree Method) -C ---------------------------------------- - CALL ORDRB5(S) - IF(.NOT.DEBUGSPK) MSGLVA=0 -C - RETURN - END ##### sed > DualAffCode/main.f <<'#####' 's/^-//' -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 * It uses, at each iteration, some subroutines of the * -C * SPARSPAK-A software package to find the direction of * -C * improvement. Also uses linpack routines... * -C *************************************************************** - INCLUDE (common.h) - CHARACTER*70 argv - DOUBLE PRECISION W(MXRWS), V(MXRWS), S(MXSNZ), ATY(MXCLS) - LOGICAL KCNTCL, KDUAL, KPREPR -C - kscale=.FALSE. - kschur=.FALSE. - kcntcl=.FALSE. - ksplit=.FALSE. - kdual =.FALSE. - kprepr=.FALSE. - IDNTHR=0 - iargcntr=1 - 2 if (iargcntr.gt.iargc()) goto 4 - call getarg(iargcntr,argv) - if( argv.eq.'-d' ) kschur=.TRUE. - if( argv.eq.'-t' ) ksplit=.TRUE. - if( argv.eq.'-s' ) kscale=.TRUE. - if( argv.eq.'-u' ) kdual=.TRUE. - if( argv.eq.'-p' ) kprepr=.TRUE. - if( argv.eq.'-c' ) kcntcl=.TRUE. - if( argv.eq.'-h' .or. argv.eq.'-?' ) then - write(0,*) ' da [-s] [-c] [-u] [-p] [-d ###] [-t ###] =### nzs' - write(0,*) ' -t dense cols have avg. rows bet. nz >=###' - stop - endif - if( argv.eq.'-d'.or.argv.eq.'-t' ) then - iargcntr=iargcntr+1 - call getarg(iargcntr, argv) - read(argv,*) IDNTHR - write(6,*) ' Using Schur Complement: idnthresh: ',idnthr - endif - iargcntr=iargcntr+1 - goto 2 - 4 continue -C -------------------------- -C Initialize SPARSPAK-A -C -------------------------- - CALL SPRSPK - MAXSA=MXSNZ -C --------------------- -C INITIALIZE PARAMETERS -C --------------------- - ITIM(1)=MCLOCK() - CALL INTLZE -C - CALL READAT - NSTFRE= NFREE - ITIM(2)=MCLOCK() - CALL SORTCL - IF (KSCALE) CALL SCALE - IF (KCNTCL.OR.KSPLIT) THEN - CALL COLCOUNT - IF(KCNTCL) STOP - ENDIF - IF(KPREPR) CALL PREPRO - IF(KDUAL ) CALL DUALIZ - ITIM(3)=MCLOCK() - IF (FEASBL .NE. 'FEASIBLE') THEN - WRITE(0,*) ' After Preproc: FEASBL: ',FEASBL - GOTO 170 - ENDIF - CALL FREEVA - CALL ISLACK - CALL STARTY(2) -C................................................... Calc. Slacks - CALL DCOPY(NSPCOLS,C,1,SLACK,1) - CALL MULATY(Y,ATY, ISUBN) - CALL DAXPY(NSPCOLS,-DBONE,ATY,1,SLACK,1) -C .............................................. Find slkmin - SLKMIN=UBSTEP - DO 10 I=1,NSPCOLS - IF (SLACK(I). LE. SLKMIN) SLKMIN = SLACK(I) - 10 CONTINUE - IF(SLKMIN . GT. -DBONE .AND. SLKMIN .LE. ZEROD) SLMIN = -DBONE - IF (SLKMIN .LE. ZEROD) THEN - PHASE=1 -C .......................................... Set Phase 1 - BIGM = BIGMIU* DABS( DDOT(NROWS,B,1,Y,1) )/(-2.0D0 * SLKMIN) - write(6,*) ' BIGM: ',BIGM - Y(NROWS+1) = (-2.0D0*SLKMIN) - B(NROWS+1) = -BIGM -C ......................................... ReCalc. Slacks - CALL DCOPY(NSPCOLS,C,1,SLACK,1) - CALL MULATY(Y,ATY,ISUBY) - CALL DAXPY(NSPCOLS,-DBONE,ATY,1,SLACK,1) - ELSE - PHASE=2 - ENDIF -C ..................... NOTE: SPLITDENSE does nothing if .NOT.KSCHUR - CALL SPLITDENSE - if(ksplit) kschur=.true. - CALL ARWISE - CALL INIMTR - CALL OUTPRD(JERROR) - IF (JERROR .GT. ZEROI) GOTO 170 - CALL INITSP(S) - ITIM(4)=MCLOCK() -C - WRITE(6,800) - DO 100 KK=1,MAXITS - ITIM(10)=ITIM(10)-MCLOCK() - CALL FINDDY(S) - ITIM(10)=ITIM(10)+MCLOCK() - ITIM(11)=ITIM(11)-MCLOCK() - CALL MULATY(DY,DV, ISUBP ) - CALL DSCAL(NCOLS,-DBONE,DV,1) - CALL STEPSZ - ITIM(11)=ITIM(11)+MCLOCK() - IF (FEASBL.NE.'FEASIBLE') GOTO 120 - CALL DAXPY(NROWS,STEP,DY,1,Y,1) - IF (PHASE .EQ. 1) THEN - Y(NROWS+1)= Y(NROWS+1) + STEP*DY(NROWS+1) - IF (Y(NROWS+1) .GE. ZEROD) THEN - IF((YOLD - Y(NROWS+1)) .GE. EPSSTP) THEN - YOLD=Y(NROWS+1) - GOTO 95 - ELSE IF(Y(NROWS+1) .GT. EPSPH1) THEN - FEASBL= 'STUCK-2' - WRITE(0,*) 'Stuck2:(delta(y))=',YOLD - Y(NROWS+1) - WRITE(0,*) ' ART: ',Y(NROWS+1) - ZLOWER=DDOT(NROWS,B,1,Y,1) - GOTO 120 - ELSE - FEASBL= 'STUCK-3' - WRITE(0,*) 'Stuck3:(delta(y))=',YOLD - Y(NROWS+1) - WRITE(0,*) ' ART: ',Y(NROWS+1) - ZLOWER=DDOT(NROWS,B,1,Y,1) - GOTO 120 - ENDIF - ELSE IF(Y(NROWS+1).LT.ZEROD) THEN - ITERA1=ITERTN - ITIM(5)=MCLOCK() - PHASE=2 - ENDIF - ENDIF - ZLOWER=DDOT(NROWS,B,1,Y,1) - IF(ZLOWER.GT.UBOUND) THEN - FEASBL='UNBOUNDD' - GOTO 120 - ENDIF -C ...................................................... Stopping Criterion 1 - AUX=DABS(ZOLD) - IF (AUX .LT. DBONE) AUX=DBONE - AUX2=DABS(ZLOWER-ZOLD)/AUX - IF (AUX2 .LE. EPSSTP) GOTO 120 - ZOLD=ZLOWER -C ...................................................... Stoppinf Criterion 2 - 95 CONTINUE - CALL GETSOL(3) -c AUX2= DABS( ZUPPER - ZLOWER ) -c IF ( AUX2 .LE. EPSSTP ) GOTO 120 -C .......................................... Update for next itertn. - IF (.NOT. KSCHUR) WRITE(6,801) KK, AUX2, Y(NROWS+1), ZLOWER, - 1 RCDNO(1) - IF (KSCHUR) WRITE(6,802) KK, AUX2, Y(NROWS+1), ZLOWER, - 1 RCDNO(1), RCDNO(2) -C ......................................... ReCalc. Slacks and updmtr - ITIM(12)=ITIM(12)-MCLOCK() - CALL DCOPY(NCOLS,C,1,SLACK,1) - CALL MULATY(Y,ATY, ISUBP) - CALL DAXPY(NCOLS,-DBONE,ATY,1,SLACK,1) - CALL UPDMTR - ITIM(12)=ITIM(12)+MCLOCK() - 100 CONTINUE -C - 120 IF (PHASE .EQ. 2) THEN - ITERA2= KK-ITERA1 - ITIM(6)=MCLOCK() - ELSE - ARTFCL= Y(NROWS+1) - ITERA1= KK - ITIM(5)=MCLOCK() - ENDIF - CALL GETSOL(-1) - CALL PRINT -C - 170 STOP - 800 FORMAT('ITER',7X,'STOPC',9X,'YART',9X,'ZLOWER',9X,'ChCond', - 1 9X,'ScCond') - 801 FORMAT(1x,I3,2X,4(E13.5,1X)) - 802 FORMAT(1x,I3,2X,5(E13.5,1X)) -C - END ##### sed > DualAffCode/makefile <<'#####' 's/^-//' -.SUFFIXES : .o $(.SUFFIXES) - -OBJECTS = bldupd.o finddy.o numerics.o read.o wrapup.o debug.o main.o prep.o tricks.o - -OPTIONS = -O - -.f.o : - xlf -c $(OPTIONS) $< - -da1 : $(OBJECTS) - xlf $(OPTIONS) -o da $(OBJECTS) /usr/um/math/lib/linpack_d.a -L/afs/engin.umich.edu/group/engin/ioe/spk/AIX -lspk - -all: da1 - -# Note: to compile on non RS/6000 machines, use f77 instead of xlf. ##### sed > DualAffCode/notes <<'#####' 's/^-//' ------------------------------------------------------------ - Dual Affine code (DA) notes: by dfholmes, 7/29/92 ------------------------------------------------------------ - -I. General notes. - - The program da is a dual affine code described in Adler, et. -al. in ORSA journal on comuting vol.1, no. 2. It is generally -designed to be easily modifiable so that other algorithms (e.g. -primal-dual) affine code. may also be easily implemented and -tested. - - The code was originally written by Jose Arantes, now of the -University of Cincinnati, and has been ported to unix and substantially -modified by D.F. Holmes. - -To run the da code, you need an mps file describing the linear -program. Type da [options] output file. To get -help, type da -h. Doing so produces: - - da [-s] [-c] [-u] [-p] [-d ###] [-t ###] =### nzs - -t dense cols have avg. rows bet. nz >=### -STOP -wheel% - -Full descriptions of the options are: - -[-s] Scale problem. da uses the same scaling algorithm that OB1 - uses, i.e. (1) divide each nonzero by row norm, (2) divide - each nonzero by its column norm, (3) repeat (1) and (2), - (4) divide each row by its largest absolute value, and (5) - divide each column by its largest absolute value. - -[-d ###] Use Schur complement. Declare all rows with no less than - ### nonzeroes as dense. - -[-c] Give nonzero report: Reads in mps file, and prints - (1) column number, (2) column name, (3) nonzeros in column, - (4) average number of rows between nonzero elements, and - (5) column density (no. of nonzeroes/no. of rows) - -[-u] Form dual of problem in internal data structures before - solving. This hasn't been fully tested since being ported to - UNIX. - -[-p] Preprocess problem as described in Adler, et. al. (see above - for referencE). This hasn't been tested since being ported to - UNIX. - -[-d ###] Use Schur complement. Declare all rows with no less than - ### nonzeroes as dense. - -[-t ###] Use Schur complement. Declare all cols with an average - number of rows between nonzeros no less than ### to be dense. - - -II. Notes on programing and modification. - - The program has been split up into the following files, each with one -or more fortran subroutines. - -main.f : Main program and interface to UNIX. -finddy.f : Contains interface code to SPARSPAK and routines to calc. - search directions. -bldupd.f : Contains routines to build the search direction matrices (ADAt) - and to update them as the algorithm progresses. -numerics.f:Contains support numerical routines. -tricks.f : Implements the options described above. -prep.f : Contains preprocessing and initialization routines. -read.f : Contains code to read in the MPS format and store it - using sparse array techniques. Parts of this code are still messy. -wrapup.f : Contains subroutines to print final solution, data, etc. -debug.f : Contains a subroutine that prints out parts of the data structure. - Use calls to this routine (prtdeb) to find out what's going - on as the algorithm progresses. - -Almost all data is passed to each subroutine through a master include file -called "common.h". To change array sizes, etc. one only need change the -parameters in common.h, and recompile everything. A makefile is provided -to ease compilation; just type "make" to compile all of the source code. - -Two user libraries are needed to compile this code; (1) SPARSPAK, from George` -and Liu, (available as a library as /afs/engin.umich.edu/group/engin/ioe/spk/libspk.a) -and LINPACK, available as /usr/um/math/lib/linpack_d.a). References to these -are already included in the makefile. To date, only the RS/6000s have sparspak, -but it should be easy to make them for other platforms. - -To implement other algorithms, you'll have to change INIMTR and UPDMTR in -bldupd.f and FINDDY in finddy.f. (maybe others). The subroutine FINDDG -in tricks.f can be modified to calculate new values for diagonal elements -in the search direction system... (Gee, isn't it easy?) - -III. Notes on schur complement processing: - -General approach: Read in data, decide what's dense, split these columns into a -separate (sparse) data structure, and use LINPACK/SPARSPAK to find the search -direction. All interfaces with long vectors (e.g. c, b, etc.) are done -using MPSP2C and MPDN2C (see below). Deciding if a column is dense is -arbitrary; just create code that sets IDFLAG (again, see below) before the -SPLITDENSE subroutine is processed. In general, I SPLIT AFTER preprocessing -and scaling. - -Other notes: - (1) Only the sparse part of ADAt is needed for the cholesky factorization, -so only sparse columns need to be updated. - - (2) If an eliminated dense column takes all of a given rows' nonzeros with -it, then the cholesky factorization will be ill-conditioned. In INIMTR, these -rows show up as zero pivots. I restore full row rank by setting these diagonal -entries as epsilon. Setting them to One worked for a few problems, but not for -others. - -Relevant/ New Data structures. - IAD,JAD,AD: Sparse represenatation of dense columns. - IADT, JATD, ADT: Sparse representation of transpose of dense columns. - IDNSTHR : Dense threshold. - MPSP2C(I) : Gives the actual column index for column I in (IA,JA,A) - MPDN2C(I) : Gives the actual column index for column I in (IAD,JAD,AD) - IDFLAG(I) : An alternative: If true, split into schur complement. - If false, don't. ##### sed > DualAffCode/numerics.f <<'#####' 's/^-//' -C *************************************************************** -C * Subroutine GETSOL: Gets solution -C * IOPT: 1 = GETSTX 2=GETZ 3=PRMLOB 0=NONE -1=ALL -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 GETSOL(IOPT) - INCLUDE (common.h) -C - DOUBLE PRECISION AUX(MXCLS) -C - if(iopt.eq.1.or.iopt.lt.0) then -C ............................................. GETSTX - 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 - K1= ZEROI - K2= ZEROI - DO 50 J=1,NCOLS - 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) - CALL DROPCL(K2) - 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 CONTINUE - endif - if(iopt.eq.2.or.iopt.lt.0) then -C .......................................................... GETZ -C - IF (PROBLM .EQ. 'DUAL') GOTO 190 - DUMMY= ZUPPER - ZUPPER= -ZLOWER - ZLOWER= -DUMMY -C - 190 ZUPPER= ZUPPER + ZCONST - ZLOWER= ZLOWER + ZCONST -C - endif - if(iopt.eq.3.or.iopt.lt.0) then -C ......................................................... PRMLOB - ZUPPER=ZEROD - DO 200 I=1,NCOLS - X(I)= -DBONE*DV(I)/(SLACK(I)*SLACK(I)) -c IF (X(I) .LT. ZEROD) THEN -c X(I)=-X(I) -c ENDIF -c ZUPPER= ZUPPER+C(I)*X(I)/COLSF(I) - ZUPPER= ZUPPER+C(I)*X(I) - 200 CONTINUE - endif -C - RETURN - END -C *********************************************************** -C * Subroutine mulaty(Y,ATY, isubnr) multiplies the (nxm) A'* -C * matrix specified in a dense storage scheme by a dense * -C * vector y (m x 1). The vec obtained is stored in ATY * -C * if ISUBNR= 1|| ISUBNR=-1 and phase=1 sub. dvec(nrows+1) * -C * modified 6/21/92 dfholmes to take care of schur comple. * -C *********************************************************** - SUBROUTINE MULATY(DVEC,ATY, ISUBNR) - INCLUDE (common.h) -C - DOUBLE PRECISION DVEC(MXRWS),ATY(MXCLS) -C ------------------------------------------------- Sparse - DO 20 I=1,NSPCOLS - MAP = MPSP2C(I) - ATY( MAP )= ZEROD - DO 10 J=IA(I),IA(I+1)-1 - ATY(MAP)= ATY(MAP)+A(J)*DVEC(JA(J)) - 10 CONTINUE - IF(ISUBNR.EQ.ISUBY.or.(ISUBNR.EQ.ISUBP .AND. PHASE.EQ.1)) - 1 ATY(MAP)=ATY(MAP) - DVEC(NROWS+1) - 20 CONTINUE -C ------------------------------------------------- Dense - DO 40 I=1,NDNCOLS - MAP = MPDN2C(I) - ATY( MAP )= ZEROD - DO 30 J=IAD(I),IAD(I+1)-1 - ATY(MAP)= ATY(MAP)+AD(J)*DVEC(JAD(J)) - 30 CONTINUE - IF(ISUBNR.EQ.ISUBY.or.(ISUBNR.EQ.ISUBP .AND. PHASE.EQ.1)) - 1 ATY(MAP)=ATY(MAP) - DVEC(NROWS+1) - 40 CONTINUE - RETURN - END -C -C ************************************************************* -C * * -C * Subroutine STARTY, finds an initial value of Y given by * -C * (Norm(C)/Norm(At*b))*b. * -C * * -C * called by: Main Program * -C * subroutines called: MULATY, NORM * -C ************************************************************* -C - SUBROUTINE STARTY(IOPT) - INCLUDE (common.h) -C - DOUBLE PRECISION ATB(MXCLS) -C - IF(IOPT.eq.1) then - CALL MULATY(B, ATB, ISUBN) - CNRM = DNRM2(NSPCOLS,C,1) - ATBNRM = DNRM2(NSPCOLS,ATB,1) - IF (ATBNRM .LE. EPSSTY ) ATBNRM= DBONE - CALL DCOPY(NROWS,B,1,Y,1) - CALL DSCAL(NROWS,CNRM/ATBNRM,Y,1) - IF (IDEB(1).GT.0) WRITE(6,900) - ELSE IF(IOPT.eq.2) then - DO 10 I=1,NROWS - Y(I)=DBONE - 10 CONTINUE - ELSE - WRITE(0,*) ' ERROR: Invalid STARTY option given ' - STOP - ENDIF - IF(IDEB(1).GT.1) WRITE(6,950) (i,Y(I),I=1,NROWS) - RETURN -C - 900 FORMAT(' Starting w/ N(c)/N(Atb) . b ') - 950 FORMAT(5(I3,':',E12.5)) - END - -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 - INCLUDE (common.h) -C - STEP= UBSTEP - DO 10 I=1,NSPCOLS+NDPCOLS - IF (DV(I) .GE. ZEROD) GOTO 10 - TEMP= -SLACK(I)/DV(I) - IF (TEMP .LT. STEP) STEP=TEMP - 10 CONTINUE - STEP= GAMMA*STEP - IF ((STEP/GAMMA) .GE. UBSTEP) FEASBL= 'STUCK-1 ' - RETURN - END -c stolen shamelessly from sparspak. this routine estimates the -C condition number of the cholesky factorization. I didn't -C call erest because eresti gives statistics that I don't want. -C dfholmes, 7/29/92 -C - SUBROUTINE ESTCOND ( RCOND, S ) 35SPK -C 36SPK -C*************************************************************** 37SPK -C 38SPK - INTEGER DIAG , FIRST , ICPAD , IERRA , IMPAD , 39SPK - 1 INVP , IPRNTE, IPRNTS, LINK , LNZ , 40SPK - 1 MAXINT, MAXSA , METHOD, MSGLVA, MXREQD, 41SPK - 1 MXUSED, NEDGES, NEQNS , NLONG , NLP1 , 42SPK - 1 NOFNZ , NOFSUB, NSHORT, NSP1 , NVARS , 43SPK - 1 NZSUB , PERM , RHS , STAGE , TEMPU , 44SPK - 1 TEMPV , XLNZ , XNZSUB 45SPK - REAL ALOCTM, ALOSTR, ANORM , DUMMY , ERRFCT, 46SPK - 1 ERROPS, ERRSTR, ERRTIM, FCTIME, FCTOPS, 47SPK - 1 MCHEPS, ORDSTR, ORDTIM, OVERHD, RATIOL, 48SPK - 1 RATIOS, RCONDA, RELEST, SLVOPS, SLVSTR, 49SPK - 1 SLVTIM, SVPAD , CLOCK 50SPK - DOUBLE PRECISION S(1) 51SPK - DOUBLE PRECISION DSIGMA, EMAX , EPS , FCTERR, 52SPK - 1 NORM , RCOND , RELERR, SCALE , 53SPK - 1 T , YMAX , YNORM 54SPK - DOUBLE PRECISION OPS , OPS1 , OPS2 , OPS3 , 55SPK - 1 OPS4 56SPK -C 57SPK -C*************************************************************** 58SPK -C 59SPK - COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 60SPK - 1 MCHEPS, CLOCK 61SPK - COMMON /SPAUSR/ MSGLVA, IERRA , MAXSA , NVARS 62SPK - COMMON /SPACON/ STAGE , MXUSED, MXREQD, NEQNS , NEDGES, 63SPK - 1 METHOD, NOFNZ , NOFSUB, ICPAD(38), 64SPK - 1 NSHORT, NLONG , NSP1 , NLP1 65SPK - COMMON /SPAMAP/ PERM , INVP , RHS , DIAG , XLNZ , 66SPK - 1 LNZ , XNZSUB, NZSUB , LINK , FIRST , 67SPK - 1 TEMPV , IMPAD(39) 68SPK - COMMON /SPADTA/ ORDTIM, ALOCTM, FCTIME, SLVTIM, ERRTIM, 69SPK - 1 FCTOPS, SLVOPS, ERROPS, ORDSTR, ALOSTR, 70SPK - 1 SLVSTR, ERRSTR, OVERHD, ANORM , RCONDA, 71SPK - 1 ERRFCT, RELEST, SVPAD(33) 72SPK -C 73SPK -C 81SPK - TEMPU = TEMPV + NEQNS 82SPK - MXREQD = TEMPU + NEQNS - 1 83SPK - IF ( MXREQD .GT. MAXSA ) THEN - WRITE(0,*) ' ERROR: STOP' - STOP - ENDIF - MXUSED = MXREQD 85SPK - ERRSTR = MXREQD 86SPK -C 87SPK -C -------------------------- 88SPK -C ESTIMATE CONDITION NUMBER. 89SPK -C -------------------------- 90SPK - SCALE = 0.50D0 91SPK - DUMMY = DTIME(0) 92SPK - CALL COND51 ( NEQNS, S(XLNZ), S(LNZ), S(XNZSUB), S(NZSUB), 93SPK - 1 S(DIAG), S(TEMPV), SCALE, EMAX, OPS1 ) 94SPK - CALL COND52 ( NEQNS, S(XLNZ), S(LNZ), S(XNZSUB), S(NZSUB), 95SPK - 1 S(DIAG), S(TEMPV), SCALE, YNORM, YMAX, EMAX, 96SPK - 1 OPS2) 97SPK - CALL COND53 ( NEQNS, S(XLNZ), S(LNZ), S(XNZSUB), S(NZSUB), 98SPK - 1 S(DIAG), S(TEMPV), SCALE, YNORM, OPS3 ) 99SPK - CALL COND54 ( NEQNS, S(XLNZ), S(LNZ), S(XNZSUB), S(NZSUB), 100SPK - 1 S(DIAG), S(TEMPV), SCALE, YNORM, OPS4 ) 101SPK - OPS = OPS1 + OPS2 + OPS3 + OPS4 102SPK -C 103SPK - RCOND = 0.0D0 104SPK - NORM = ANORM 105SPK - IF ( NORM .EQ. 0.0D0 ) GO TO 100 106SPK - RCOND = DMIN1(YNORM/NORM,EMAX/(NORM*YMAX)) 107SPK - OPS = OPS + 3.0D0 108SPK -C 109SPK - 100 CONTINUE 110SPK - T = 1.0D0 + RCOND 111SPK - IF ( T .NE. 1.0D0 ) GO TO 200 112SPK - RCOND = -1.0D0 113SPK -C 115SPK - 200 CONTINUE 116SPK - RETURN - END - ##### sed > DualAffCode/prep.f <<'#####' 's/^-//' -C ************************************************************ -C * Subroutine intlze sets the values of the parameters * -C * and the initial values of some variables. * -C ************************************************************ - SUBROUTINE INTLZE - INCLUDE (common.h) -C - BIGMIU= 1.0D+05 - BLANK1= ' ' - BLANK2= ' ' - DBONE = 1.0D+00 - EPSILO= 1.0D-06 - EPSSTP= 1.0D-06 - EPSPH1= 1.0D-06 - EPSUPD= 0.1D0 - EPSSTY= 1.0D-06 - ZEROI = 0 - ZEROD = 0.0D0 - FEASBL= 'FEASIBLE' - GAMMA = 0.95D0 - IALG = IDUALA - ITERA1= ZEROI - ITERA2= ZEROI - MAXITS= 100 - NFIX = ZEROI - NFREE = ZEROI - NCHSGN= ZEROI - NLBD = ZEROI - NRWDRP= ZEROI - NSTVAR= ZEROI - NDNCOLS= ZEROI - PHASE = 2 - PROBLM= 'DUAL' - SLKMIN= 1.0D+05 - UBSTEP= 1.0D+20 - UBOUND= 1.0D+20 - YOLD = UBOUND - ZCONST= ZEROD - ZOLD = -1.0D+20 - DO 10 I=1,30 - ITIM(I)=0 - 10 CONTINUE -C - RETURN - END -C -C ******************************************************************** -C * Subroutine Chkfre, verifies if a given variable is or not unres-* -C * tricted in sign. If yes, the parameter FRE returns the position* -C * occuppied by this variable in the array VRFREE. * -C * Chkfre checks all columns, even if dense. * -C ******************************************************************** - FUNCTION ICHKFRE(ICOLUMN) - INCLUDE (common.h) -C - ICHKFRE= ZEROI - DO 10 I=1,NFREE - IF (ICOLUMN .EQ. VRFREE(I)) THEN - ICHKFRE= I - RETURN - END IF - 10 CONTINUE - RETURN - END -C ******************************************************************** -C * Subroutine Chkrws, checks if there are still variables in a cons* -C * traint after a variable is fixed. If not, checks if the row is * -C * consistant and drop it if yes, otherwise the problem is infeas. * -C ******************************************************************** - SUBROUTINE CHKRWS(LESSRW) - INCLUDE (common.h) -C - INTEGER ISPNZ -C - I=1 - 10 ISPNZ= IAT(I+1)-IAT(I) - IF (ISPNZ .GT. ZEROI) GOTO 30 - IF (B(I) .GT. ZEROD) THEN - IF (RTYPE(I+1) .EQ. ' L ') GOTO 20 - FEASBL= NAMERW(I+1) - RETURN - ELSE - IF (B(I) .EQ. ZEROD) GOTO 20 - IF (RTYPE(I+1) .EQ. ' G ') GOTO 20 - FEASBL= NAMERW(I+1) - RETURN - ENDIF -C - 20 CALL DROPRW(I) - LESSRW= LESSRW+1 - CALL UPTRWD(I) - GOTO 40 -C - 30 I=I+1 - 40 IF (I .LE. NROWS) GOTO 10 -C - RETURN - END -C ******************************************************************** -C * Subroutine Chksgn, the parameter FRE returns the value zero if * -C * the given constraint does NOT have any free variable and all the * -C * nonzero coefficients have the same sign. Otherwise the returned * -C * value of FRE is POSITIVE. dfh: not checked yet. * -C ******************************************************************** - FUNCTION ICHKSGN(IROW) - INCLUDE (common.h) -C ................................................. Check Sparse ..... - ICHKSGN=ZEROI - DO 10 I= IAT(IROW), IAT(IROW+1) - 2 - IF (AT(I)*AT(I+1) .LT. ZEROD) THEN - ICHKSGN=9999 - RETURN - ENDIF - DTEMP=AT(I)*AT(I+1) - 10 CONTINUE - DO 20 I= IAT(IROW), IAT(IROW+1)-1 - ICHKSGN=ICHKFRE(ORGIDX(JAT(I))) - IF (ICHKSGN .GT. ZEROI) RETURN - 20 CONTINUE - 40 RETURN - END -C **************************************************************** -C * Subroutine Prepro: Preprocessing the Input Data. Basically, * -C * (1) Remove all rows that have all the coefficient entries * -C * with the same sign and RHS=0, set all variables that * -C * appear in these constraints to zero. * -C * (2) Declare the problem infeasible if a row is found such * -C * that all coefficient entries have the same sign and * -C * rhs has the opposite sign. * -C * (3) Rows with only one entry, fix the value of the corres- * -C * ponding variable, or change the lower bound or declare * -C * the problem infeasible, depending on the row sign. * -C **************************************************************** -C - SUBROUTINE PREPRO - INCLUDE (common.h) -C - INTEGER ROW, CLBNDD, TOTNZR, FRE, COL -C - 10 I=1 - LESSCL= ZEROI - LESSRW= ZEROI -C - 20 IF (IAT(I) .NE. (IAT(I+1)-1)) GOTO 140 - AUX= B(I)/AT(IAT(I)) - COL= JAT(IAT(I)) - ORGCOL= ORGIDX(COL) - IF (RTYPE(I+1) .NE. ' E ') GOTO 50 - FRE=ICHKFRE(ORGCOL) - IF (FRE .EQ. ZEROI) THEN - IF (AUX .LT. ZEROD) GOTO 150 - GOTO 30 - END IF -C ------------------------------------------------- -C The current variable is unrestricted in sign. -C ------------------------------------------------- - CALL UPTVFR(FRE) -C - 30 CALL FIXVAR(AUX,COL) - LESSCL= LESSCL+1 - CALL UPTCLD(COL) - CALL CHKRWS(LESSRW) - GOTO 210 -C - 40 ROW=I - CALL DROPRW(ROW) - LESSRW= LESSRW+1 - CALL UPTRWD(ROW) - GOTO 210 -C - 50 IF (IA(COL) .NE. (IA(COL+1)-1)) GOTO 60 - CALL ONEONE(AUX,I,COL,ORGCOL) - IF (FEASBL .EQ. 'FEASIBLE') GOTO 30 - GOTO 220 -C - 60 IF (AT(IAT(I)) .LT. ZEROD) GOTO 100 - FRE=ICHKFRE(ORGCOL) - IF (FRE .NE. ZEROI) GOTO 70 - IF (B(I) .GT. ZEROD) GOTO 200 - IF (B(I) .LT. ZEROD) GOTO 150 - AUX= ZEROD - GOTO 30 -C - 70 CALL UPTVFR(FRE) - CALL CHSIGN(AUX,NAMECL(COL)) -C -------------------------------- -C changing sign in the array AT -C -------------------------------- - DO 90 I1= IA(COL), IA(COL+1)-1 - ROW= JA(I1) - DO 80 I2= IAT(ROW), IAT(ROW+1)-1 - IF (JAT(I2) .EQ. COL) THEN - AT(I2) = -AT(I2) - GOTO 90 - END IF - 80 CONTINUE - 90 CONTINUE - GOTO 120 -C - 100 FRE=ICHKFRE(ORGCOL) - IF (FRE .NE. ZEROI) GOTO 110 - IF (B(I) .LT. ZEROD) GOTO 130 - GOTO 40 -C - 110 CALL UPTVFR(FRE) - 120 IF (B(I) .EQ. ZEROD) GOTO 40 - 130 CALL LOBOVA(AUX,NAMECL(COL)) - GOTO 40 -C - 140 FRE= ICHKSGN(I) - IF (FRE .GT. ZEROI) GOTO 200 - IF (B(I) .GE. ZEROD) GOTO 160 - IF (AT(IAT(I)) .LT. ZEROD) GOTO 200 - 150 FEASBL= NAMERW(I+1) - GOTO 220 -C - 160 IF (AT(IAT(I)) .LT. ZEROD) GOTO 170 - IF (B(I) .EQ. ZEROD) GOTO 180 - GOTO 200 -C - 170 IF (RTYPE(I+1) .EQ. ' L ') THEN - CALL ONENZR(I,IAT(I+1)-1,LESSCL) - GOTO 40 - END IF - IF (B(I) .NE. ZEROD) GOTO 150 -C - 180 TOTNZR=IAT(I+1) - IAT(I) - DO 190 I1=1,TOTNZR - COL= JAT(IAT(I)) - CALL FIXVAR(ZEROD,NAMECL(COL)) - LESSCL= LESSCL+1 - CALL UPTCLD(COL) - 190 CONTINUE - CALL CHKRWS(LESSRW) - GOTO 210 -C - 200 I=I+1 - 210 IF (I .LE. NROWS) GOTO 20 - IF ((LESSCL .GT. ZEROI) .OR. (LESSRW .GT. ZEROI)) GOTO 10 -C - 220 RETURN - END -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 * Remark: each row in JAT is sorted in crescent order of * -C * column index ( a result of the ARWISE subrout.)* -C ********************************************************** -C - SUBROUTINE ONENZR(IROW,ENDROW,LESSCL) - INCLUDE (common.h) -C - INTEGER BEGROW, ENDROW, BEGCOL, ENDCOL, COL , LESSCL -C - BEGROW=IAT(I) - NCOLFX=ZEROI - DO 10 I=BEGROW,ENDROW - COL= JAT(I) - NCOLFX - BEGCOL= IA(COL) - ENDCOL= IA(COL+1)-1 - IF (BEGCOL .NE. ENDCOL) GOTO 10 - IF (C(COL) .LT. ZEROD) THEN - FEASBL= NAMECL(COL) - GOTO 20 - END IF - CALL FIXVAR(ZEROD,NAMECL(COL)) - NCOLFX= NCOLFX+1 - 10 CONTINUE -C - LESSCL= LESSCL+NCOLFX -C - 20 RETURN - END -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 - SUBROUTINE ONEONE(AUX,IPTR,COL,ORGCOL) - INCLUDE (common.h) -C - INTEGER COL, ORGCOL, ROW, FRE -C - ROW=I - COEF=AT(IAT(I)) - COST=C(COL) - IF (COEF .GT. ZEROD) GOTO 10 - IF (COST .LT. ZEROD) GOTO 30 - FRE=ICHKFRE(ORGCOL) - IF (FRE .NE. ZEROI) THEN - CALL UPTVFR(FRE) - ELSE - IF (AUX .LT. ZEROD) THEN - AUX= ZEROD - END IF - END IF - GOTO 40 -C - 10 FRE=ICHKFRE(ORGCOL) - IF (FRE .EQ. ZEROI) GOTO 20 - IF (COST .GT. ZEROD) GOTO 30 - CALL UPTVFR(FRE) - GOTO 40 -C - 20 IF (AUX .LT. ZEROD) THEN - FEASBL= NAMERW(ROW+1) - GOTO 40 - END IF - IF (COST .GT. ZEROD) THEN - AUX= ZEROD - END IF - GOTO 40 -C - 30 FEASBL= NAMECL(COL) -C - 40 RETURN - END -C -C -C ********************************************************* -C * Subroutine Dropcl, eliminates a column of the LP by * -C * updating NSPCOLS and the arrays C, NAMECL, * -C * IA, JA, A AND ORGIDX. * -C ********************************************************* -C - SUBROUTINE DROPCL(ICOL) - INCLUDE (common.h) -C - INTEGER TOTNZC -C - TOTNZC= IA(ICOL+1)-IA(ICOL) - NSPCOLS= NSPCOLS-1 -C - IF (ICOL .EQ. (NSPCOLS+1)) GOTO 25 - DO 10 I=ICOL,NSPCOLS - 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 - DO 20 I= IA(ICOL), IA(NSPCOLS+1)-1 - JA(I)= JA(I+TOTNZC) - A(I)= A(I+TOTNZC) - 20 CONTINUE -C - 25 DO 30 I= IA(NSPCOLS+1), IA(NSPCOLS+2) - JA(I)= ZEROI - A(I)= ZEROD - 30 CONTINUE -C - IA(NSPCOLS+2)= ZEROI - C(NSPCOLS+1)= ZEROD - COLSF(NSPCOLS+1)=DBONE - NAMECL(NSPCOLS+1)=BLANK2 -C - RETURN - END -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 ********************************************************* - SUBROUTINE DROPRW(ROW) - INCLUDE (common.h) -C - INTEGER ROW, RWDROP -C - NRWDRP= NRWDRP+1 - RWDROP=ROW - NROWS= NROWS-1 - DO 10 I= RWDROP,NROWS - RTYPE(I+1)= RTYPE(I+2) - NAMERW(I+1)= NAMERW(I+2) - B(I)= B(I+1) - ROWSF(I)=ROWSF(I+1) - 10 CONTINUE - RTYPE(NROWS+2)= BLANK1 - NAMERW(NROWS+2)= BLANK2 - B(NROWS+1)= ZEROD - ROWSF(NROWS+1) = DBONE -C - NEXTJA= ZEROI - DO 50 I=1,NSPCOLS - IA(I)= NEXTJA+1 - DO 40 J= IA(I), IA(I+1)-1 - IF (JA(J) .EQ. RWDROP) GOTO 40 - IF (JA(J) .LT. RWDROP) GOTO 30 - JA(J)= JA(J)-1 - 30 NEXTJA= NEXTJA+1 - JA(NEXTJA)= JA(J) - A(NEXTJA)= A(J) - 40 CONTINUE - 50 CONTINUE -C - DO 60 I= NEXTJA+1, IA(NSPCOLS+1) - JA(I)= ZEROI - A(I)= ZEROD - 60 CONTINUE -C - IA(NSPCOLS+1)= NEXTJA+1 -C - RETURN - END -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 ********************************************************* - SUBROUTINE CHSIGN(VALUE,CLNAME) - INCLUDE (common.h) -C - CHARACTER*8 CLNAME -C - NCHSGN= NCHSGN+1 - IDX=INDCOL(CLNAME) - VRCSGN(NCHSGN)= ORGIDX(IDX) - DO 10 I= IA(IDX), IA(IDX+1)-1 - A(I)= -DBONE*A(I) - 10 CONTINUE - C(IDX)= -DBONE*C(IDX) - VALUE= -DBONE*VALUE -C - RETURN - END - -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 ********************************************************* - SUBROUTINE FIXVAR(VALUE,ICOL) - INCLUDE (common.h) -C - CHARACTER*8 CLNAME -C - CLNAME=NAMECL(ICOL) - IDX=INDCOL(CLNAME) - NFIX= NFIX+1 - VARFIX(NFIX)= ORGIDX(IDX) - VALFIX(NFIX)= VALUE -C ---------------- -C UPDATE RHS -C ------------------ - DO 10 I=IA(IDX), IA(IDX+1)-1 - B(JA(I))= B(JA(I)) - A(I)*VALUE - 10 CONTINUE -C --------------------------------------------- -C Update the constant of the objective function -C ---------------------------------------------- - ZCONST= ZCONST + C(IDX)*VALUE - CALL DROPCL(IDX) - RETURN - END -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 ********************************************************** - SUBROUTINE INSERW(TYPE,NAME,NZCOLS,NZVALS,RHS,TOTNZR) - INCLUDE (common.h) -C - CHARACTER*3 TYPE - CHARACTER*8 NAME -C - DOUBLE PRECISION RHS, AUXD, NZVALS(MXCLS) - INTEGER NZCOLS(MXCLS), TOTNZR -C ---------------------------------------------------------- -C IAOLD: points to the last position in JA ( or A) not yet -C updated. -C NEXTJA: points to the next position in the new arrays JA -C (or A) which will be updated. -C ----------------------------------------------------------- - NROWS= NROWS+1 - RTYPE(NROWS+1)= TYPE - NAMERW(NROWS+1)= NAME - B(NROWS)= RHS - ROWSF(NROWS)= DBONE -C - NEXTCL= NSPCOLS+1 - IAOLD= IA(NSPCOLS+1) - 1 - IA(NSPCOLS+1)= IA(NSPCOLS+1) + TOTNZR - NEXTJA= IA(NEXTCL) -1 -C ------------------------------------------------------------- -C Sorting NZCOLS and NZVALS in decreasing order of column index. -C -------------------------------------------------------------- - IF (TOTNZR .EQ. 1) GOTO 30 - DO 20 I=1,TOTNZR-1 - DO 10 J= I+1,TOTNZR - IF (NZCOLS(I) .GT. NZCOLS(J)) GOTO 10 - IF (NZCOLS(I) .EQ. NZCOLS(J)) GOTO 100 - IAUXI = NZCOLS(I) - AUXD = NZVALS(I) - NZCOLS(I)= NZCOLS(J) - NZVALS(I)= NZVALS(J) - NZCOLS(J)= IAUXI - NZVALS(J)= AUXD - 10 CONTINUE - 20 CONTINUE -C - 30 DO 90 I= 1,TOTNZR - NEXTCL= NEXTCL-1 - DO 50 J= NEXTCL,1,-1 - IF (NZCOLS(I) .EQ. J) GOTO 70 - DO 40 K= IAOLD,IA(J),-1 - JA(NEXTJA)= JA(K) - A(NEXTJA)= A(K) - NEXTJA= NEXTJA-1 - 40 CONTINUE - IAOLD= IA(J) - 1 - IA(J)= NEXTJA+1 - 50 CONTINUE - WRITE(6,60) - 60 FORMAT('ERROR. Not all the NZR coefs. were entered.') -c - 70 JA(NEXTJA)= NROWS - A(NEXTJA)= NZVALS(I) - NEXTJA= NEXTJA-1 - NEXTCL= J - DO 80 K= IAOLD,IA(J),-1 - JA(NEXTJA)= JA(K) - A(NEXTJA)= A(K) - NEXTJA= NEXTJA-1 - 80 CONTINUE - IAOLD= IA(NEXTCL)-1 - IA(NEXTCL)= NEXTJA+1 - 90 CONTINUE - GOTO 120 -C - 100 WRITE(6,110) - 110 FORMAT('ERROR. Two coefs. for the same variable.') -C - 120 RETURN - END -C ********************************************************* -C * Subroutine Lobova: it changes the lower bounded * -C * variable x >= LB to x'= x-LB >= 0. * -C * It also updates the original problem by * -C * adding a constant to the new objective * -C * function and subtracting A(I)*LB from * -C * the corresponding RHS. * -C ********************************************************* - SUBROUTINE LOBOVA(LOWBND,CLNAME) - INCLUDE (common.h) -C - CHARACTER*8 CLNAME - DOUBLE PRECISION LOWBND -C - IDX=INDCOL(CLNAME) - ZCONST= ZCONST+ C(IDX)*LOWBND - NLBD= NLBD+1 - VARLBD(NLBD)= ORGIDX(IDX) - VALLBD(NLBD)= LOWBND -C ------------------------------------ -C Modification of the Right-Hand Side -C ------------------------------------ - DO 10 I= IA(IDX), IA(IDX+1)-1 - B(JA(I))= B(JA(I)) - A(I)*LOWBND - 10 CONTINUE -C - RETURN - END -C ********************************************************** -C * Subroutine Loupva: it deals with a variable which has * -C * upper and lower bound. First it changes the * -C * variable so that the lower bound becomes * -C * zero, then it adds a nex constraint to the * -C * system in order to satisfy the upper bound. * -C ********************************************************** - SUBROUTINE LOUPVA(UPPBND,LOWBND,CLNAME) - INCLUDE (common.h) -C - CHARACTER*3 TYPE - CHARACTER*8 CLNAME, NAME - DOUBLE PRECISION UPPBND, LOWBND, NZVALS(MXCLS), RHS - INTEGER NZCOLS(MXCLS), TOTNZR -C - IF (LOWBND .EQ. ZEROD) GOTO 10 - CALL LOBOVA(LOWBND,CLNAME) -C - 10 TOTNZR=1 - IDX=INDCOL(CLNAME) - NAME= 'UPPBOUND' - NZCOLS(TOTNZR)= IDX - NZVALS(TOTNZR)= 1.0 - RHS= UPPBND - LOWBND - TYPE= ' L ' -C - CALL INSERW(TYPE,NAME,NZCOLS,NZVALS,RHS,TOTNZR) -C - RETURN - END -C ********************************************************** -C * Subroutine Freeva, finds free variables and replace * -C * them by the difference of two non-negative * -C * variables. * -C * called by: Main Program NOTENOTE: INCLUDES FREVA1 now* -C ********************************************************** - SUBROUTINE FREEVA - INCLUDE (common.h) -C - IF (PROBLM .EQ. 'DUAL') GOTO 20 - DO 10 I=1,NSPCOLS - IF (RTYPE(I+1) .NE. ' E ') GOTO 10 - IA(NSPCOLS+2)= IA(NSPCOLS+1)+IA(I+1)-IA(I) - K=IA(NSPCOLS+1)-1 - DO 8 J=IA(I),IA(I+1)-1 - K=K+1 - JA(K)=JA(J) - A(K)=-A(J) - 8 CONTINUE - NSPCOLS=NSPCOLS+1 - C(NSPCOLS)= -C(I) - COLSF(NSPCOLS)=DBONE - 10 CONTINUE - GOTO 70 -C - 20 DO 60 I=1,NFREE - DO 30 I1=1,NSPCOLS - IF (VRFREE(I) .EQ. ORGIDX(I1)) GOTO 50 - 30 CONTINUE - WRITE(6,40) VRFREE(I) - 40 FORMAT('The free variable',I4,' was not found in ORGIDX(J)') - GOTO 70 -C - 50 CONTINUE - IA(NSPCOLS+2)= IA(NSPCOLS+1)+IA(I1+1)-IA(I1) - K=IA(NSPCOLS+1)-1 - DO 58 J=IA(I1),IA(I1+1)-1 - K=K+1 - JA(K)=JA(J) - A(K)=-A(J) - 58 CONTINUE - NSPCOLS=NSPCOLS+1 - C(NSPCOLS)= -C(I1) - COLSF(NSPCOLS)=DBONE - ORGIDX(NSPCOLS)= ORGIDX(I1) - 60 CONTINUE -C - 70 RETURN - END -C ********************************************************* -C * Subroutine Islack, inserts a slack variable into an * -C * inequation to modify it to an equation. * -C ********************************************************* - SUBROUTINE ISLACK - INCLUDE (common.h) -C - DOUBLE PRECISION COEFSK -C - IF (PROBLM .EQ. 'DUAL') GOTO 30 - COEFSK= -DBONE - DO 20 I=1,NROWS - DO 10 J=1,NFREE - IF (ORGIDX(I) .EQ. VRFREE(J)) GOTO 20 - 10 CONTINUE - CALL UPTARR(COEFSK,I) - 20 CONTINUE - GOTO 100 -C - 30 DO 70 I=1,NROWS - COEFSK=ZEROD - IF (RTYPE(I+1) .EQ. ' G ') COEFSK=-DBONE - IF (RTYPE(I+1) .EQ. ' L ') COEFSK=DBONE - IF(DABS(COEFSK).GE.EPSILO) THEN - CALL UPTARR(COEFSK,I) - NAMECL(NSPCOLS)= 'SLACK ' - MPSP2C(NSPCOLS)=NSPCOLS - RTYPE(I+1)= ' E ' - ENDIF -C - 70 CONTINUE - GOTO 100 -C - 80 WRITE(6,90) - 90 FORMAT('WARNING. Adding SLACK to the objective function?') -C - 100 RETURN - END ##### sed > DualAffCode/read.f <<'#####' 's/^-//' -C *********************************************************** -C * Program Readat: reads in the data of a LP program ** -C * Written by:Jose' Arantes Jan/89 Subs. modification dfh * -C *********************************************************** - SUBROUTINE READAT - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN -C - I= ZEROI - 10 NEXTLN= BLANK1 - READ(5,'(A3)') NEXTLN - IF (NEXTLN .NE. 'ROW') THEN - I= I+1 - IF (I .EQ. 50) THEN - WRITE(6,40) - RETURN - ENDIF - GOTO 10 - ENDIF - 20 I= ZEROI - CALL READRW(NEXTLN) - IF (NEXTLN .NE. 'COL') THEN - WRITE(6,60) - RETURN - ENDIF - CALL READCL(NEXTLN) - IF (NEXTLN .NE.'RHS') THEN - WRITE(6,80) - RETURN - ENDIF - CALL READRH(NEXTLN) - IF (NEXTLN .EQ. 'RAN') CALL READRG(NEXTLN) - IF (NEXTLN .EQ. 'BOU') CALL READBD(NEXTLN) - IF (NEXTLN .EQ. 'END') RETURN - WRITE(6,100) - 100 FORMAT('ERROR: The indicator card ENDDATA was not found.') - RETURN -C - 40 FORMAT('ERROR: The indicator card ROWS was not found.') - 60 FORMAT('ERROR: The indicator card COLUMNS was not found.') - 80 FORMAT('ERROR: The indicator card RHS was not found.') - END -C********************************************************************* -C Subroutine Readrw: it reads in the row name and the (in) equality * -C sign associated to each row. The signs G, L,E correspond to cons. * -C and N to the objective function. * -C note: The obj. funct. will always be the 1st row in NAMERW,RTYPE * -C********************************************************************* - SUBROUTINE READRW(NEXTLN) - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN, AUX1 - CHARACTER*8 AUX2 -C - I= 1 - 10 RTYPE(I)= BLANK1 - NAMERW(I)= BLANK2 - READ(5,20) RTYPE(I), NAMERW(I) - 20 FORMAT(A3,1X,A8) - IF (RTYPE(I) .EQ. 'COL') GOTO 30 - I= I+1 - GOTO 10 -C - 30 NROWS= I-1 - NEXTLN= 'COL' - IF (RTYPE(1) .NE. ' N ') THEN - IOBJRW=ZEROI - DO 40 J=2,NROWS - IF (RTYPE(J) .EQ. ' N ') IOBJRW=J - 40 CONTINUE - IF (IOBJRW.EQ. ZEROI) THEN - WRITE(6,*) ' ERROR: NO OBJECTIVE ROW FOUND' - STOP - ENDIF - AUX1= RTYPE(IOBJRW) - AUX2= NAMERW(IOBJRW) - NAMERW(IOBJRW)= NAMERW(1) - RTYPE(IOBJRW)= RTYPE(1) - NAMERW(1)= AUX2 - RTYPE(1)= AUX1 - ENDIF - NROWS= NROWS-1 -C - RETURN - END -C -C********************************************************************* -C Subroutine Readcl: it reads in the cost coefficients and the non- * -C zero elements of each column, associated row index and value. It * -C also fills in the arrays C(I), IA(I), JA(I) and A(I), where: * -C C(I) has the cost coefficient of the ith variable * -C IA(I) is the position in JA and A where starts the Ith column. * -C JA(K) gives the row index of the (K-IA(J)+1)th non-zero (J=COL) * -C A(K) gives the value of the (K-IA(J)+1)th non-zero element in colJ.* -C * -C Addition (dfholmes, 6/20): Now, if column is considered dense, put * -C the column's data in (CD,IAD,JAD,AD) and set INDDNS. - LATER * -C********************************************************************* -C - SUBROUTINE READCL(NEXTLN) - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN - CHARACTER*8 CURRCL, LASTCL, ROW(MXFLDS) - DOUBLE PRECISION VAL(MXFLDS) -C - DO 5 I=1,MXCLS - C(I)= ZEROD - NAMECL(I)= BLANK2 - 5 CONTINUE - LASTCL= 'STARTCOL' - INZCTR =ZEROI - NSPCOLS=ZEROI - NDNCOLS=ZEROI -C - 10 CURRCL= BLANK2 - NEXTLN= BLANK1 - DO 15 IFIELD=1,MXFLDS - ROW(IFIELD)=BLANK2 - VAL(IFIELD)=ZEROD - 15 CONTINUE -C everything starts out as sparse - READ(5,900) NEXTLN,CURRCL,ROW(1),VAL(1),ROW(2), VAL(2) - IF (NEXTLN .EQ. 'RHS' .OR. NEXTLN .EQ. 'BOUNDS') GOTO 70 - IF (CURRCL .NE. LASTCL) THEN - NSPCOLS=NSPCOLS+1 - ORGIDX(NSPCOLS)=NSPCOLS - MPSP2C(NSPCOLS)=NSPCOLS - IA(NSPCOLS)=INZCTR+1 - LASTCL= CURRCL - NAMECL(NSPCOLS)= CURRCL - ENDIF - DO 60 IFIELD=1,MXFLDS - IF(ROW(IFIELD) .EQ. BLANK2) GOTO 60 - IF(VAL(IFIELD) .EQ. ZEROD) GOTO 60 - INDR=INDROW(ROW(IFIELD)) - IF (INDR .EQ. ZEROI) THEN - C(NSPCOLS)=VAL(IFIELD) - ELSE IF (INDR. NE. -1) THEN - INZCTR=INZCTR+1 - JA(INZCTR)= INDR - A(INZCTR) = VAL(IFIELD) - ENDIF - 60 CONTINUE - GOTO 10 -C - 70 IA(NSPCOLS+1)= INZCTR+1 - NSTVAR= NSPCOLS -C - 900 FORMAT(A3,1X,A8,2X,A8,2X,F12.4,3X,A8,2X,F12.4) - RETURN - END -C -C********************************************************************* -C Subroutine Readrh: reads in RHS * -C********************************************************************* -C - SUBROUTINE READRH(NEXTLN) - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN - CHARACTER*8 ROW(MXFLDS),CURRHS - DOUBLE PRECISION VAL(MXFLDS) -C - DO 10 I=1,NROWS - B(I)= ZEROD - 10 CONTINUE - 20 DO 25 IFIELD=1,MXFLDS - ROW(IFIELD)=BLANK2 - VAL(IFIELD)=ZEROD - 25 CONTINUE - NEXTLN= BLANK1 -C - READ(5,900) NEXTLN, CURRHS, ROW(1), VAL(1), ROW(2), VAL(2) - IF (NEXTLN .EQ. BLANK1) THEN - DO 30 IFIELD=1,MXFLDS - IF(ROW(IFIELD) .EQ. BLANK2) GOTO 30 - INDEX=INDROW(ROW(IFIELD)) - IF(INDEX.NE.-1) B(INDEX)=VAL(IFIELD) - 30 CONTINUE - GOTO 20 - ENDIF - RETURN - 900 FORMAT(A3,1X,A8,2X,A8,2X,F12.4,3X,A8,2X,F12.4) - END -C -C********************************************************************* -C Subroutine Readrg: reads in ranges for some constraint. * -C********************************************************************* -C - SUBROUTINE READRG(NEXTLN) - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN, TYPE - CHARACTER*8 ROW(MXFLDS), CURRRG, LASTRG - CHARACTER*8 NAME -C - DOUBLE PRECISION NZVALS(MXCLS), RG(MXFLDS) - INTEGER NZCOLS(MXCLS) -C - 5 CURRRG= BLANK2 - DO 7 IFIELD=1,MXFLDS - ROW(IFIELD)=BLANK2 - RG(IFIELD) =ZEROD - 7 CONTINUE - NEXTLN= BLANK1 -C - READ(5,900) NEXTLN, CURRRG, ROW(1), RG(1), ROW(2), RG(2) -C - IF (NEXTLN .NE. BLANK1) RETURN - IF (LASTRG .EQ. 'STARTRNG') THEN - WRITE(6,*) ' WARNING: RANGE REPEATED: CURRG',CURRG - GOTO 5 - ENDIF - IF (CURRRG .NE. LASTRG) THEN - LASTRG= BLANK2 - LASTRG= CURRRG - ENDIF - 30 DO 60 IFIELD=1,MXFLDS - IF(ROW(IFIELD) .EQ. BLANK2) GOTO 60 - INDR=INDROW(ROW1) -C ----------------------------------------------------------- -C Find the COLUMN indices and coefficients of non-zero entries -C of the constraint ROW1 and store them respectively in the -C arrays NZCOLS(I) and NZVALS(I), whose dimension is NCOLS. -C ------------------------------------------------------------ - DO 35 I=1,NSPCOLS - NZCOLS(I)= ZEROI - NZVALS(I)= ZEROD - 35 CONTINUE -C - KS=ZEROI - DO 50 I= 1,NSPCOLS - DO 40 J= IA(I), IA(I+1)-1 - IF (JA(J) .EQ. INDR) THEN - KS=KS+1 - NZCOLS(KS)=I - NZVALS(KS)=A(J) - ENDIF - 40 CONTINUE - 50 CONTINUE - NZSP = KS -C -------------------------------------------------------- -C NAMERW(I), RTYPE(I) and B(I) are obtained for the new -C constraint and RTYPE is updated for the old constraint. -C ------------------------------------------------------------ - RHS= B(INDR)+RG(IFIELD) - NAME= NAMERW(INDR+1) -C <= - IF (RHS .LE. B(INDR)) THEN - TYPE=' G ' - IF(RTYPE(INDR+1) .EQ. ' E ') THEN - RTYPE(INDR+1)=' L ' - ELSE IF ( RTYPE(INDR+1) .NE. ' L ' ) THEN - WRITE(6,*) ' ERROR: RANGE ',CURRRG,' INCONSISTENT' - RETURN - ENDIF - ELSE -C > - TYPE=' L ' - IF(RTYPE(INDR+1) .EQ. ' E ') RTYPE(INDR+1)=' G ' - IF(RTYPE(INDR+1) .NE. ' G ') THEN - WRITE(6,*) ' ERROR: RANGE ',CURRRG,' INCONSISTENT' - RETURN - ENDIF - ENDIF - CALL INSERW(TYPE, NAME, RHS, NZCOLS, NZVALS, NZSP) - - 60 CONTINUE - GOTO 5 - 900 FORMAT(A3,1X,A8,2X,A8,2X,F12.4,3X,A8,2X,F12.4) - END -C -C********************************************************************* -C Subroutine Readbd: it reads in the bounds restrictions on each vari* -C able and when appropriate it eliminates a variable, changes the * -C variable range and/or create a new constraint. * -C * -C Note: dfholmes. This should probably be changed in the future. * -C********************************************************************* -C - SUBROUTINE READBD(NEXTLN) - INCLUDE (common.h) -C - CHARACTER*3 NEXTLN, TYPEBO(MXCLS) - CHARACTER*8 CURRBO, LASTBO, COLBOU(MXCLS) -C - DOUBLE PRECISION VALBOU(MXCLS) -C - DO 10 I=1,MXCLS - TYPEBO(I)= BLANK1 - COLBOU(I)= BLANK2 - VALBOU(I)= ZEROD - 10 CONTINUE -C - I=1 - LASTBO= 'STARTBOU' - 20 CURRBO= BLANK2 - NEXTLN= BLANK1 -C - READ(5,30) NEXTLN, CURRBO, COLBOU(I), VALBOU(I) - 30 FORMAT(A3,1X,A8,2X,A8,2X,F12.4) -C - IF (NEXTLN .EQ. 'END') GOTO 70 - IF (CURRBO .EQ. LASTBO) GOTO 60 - IF (LASTBO .EQ. 'STARTBOU') GOTO 50 - WRITE(6,40) - 40 FORMAT('WARNING. It is provided more than one BOUND set.') - GOTO 20 -C - 50 LASTBO= BLANK2 - LASTBO= CURRBO - 60 TYPEBO(I)= NEXTLN - I= I+1 - GOTO 20 -C - 70 ITOTBO= I-1 -C - DO 300 I=1,ITOTBO - IF (TYPEBO(I) .EQ. BLANK1) GOTO 300 - IF (TYPEBO(I) .EQ. ' FX') GOTO 170 - IF (TYPEBO(I) .EQ. ' FR') GOTO 190 - IF (TYPEBO(I) .EQ. ' PL') GOTO 290 - DO 80 J=I+1,ITOTBO - IF (COLBOU(I) .EQ. COLBOU(J)) GOTO 110 - 80 CONTINUE - IF (TYPEBO(I) .EQ. ' LO') GOTO 210 - IF (TYPEBO(I) .EQ. ' UP') GOTO 220 - IF (TYPEBO(I) .EQ. ' MI') GOTO 190 -C - 90 WRITE(6,100) TYPEBO(I) - 100 FORMAT('ERROR. The value of TYPEBO(I) is ',A4) - GOTO 310 -C - 110 IF (TYPEBO(I) .EQ. ' UP') GOTO 120 - IF (TYPEBO(I) .EQ. ' LO') GOTO 130 - IF (TYPEBO(I) .EQ. ' MI') GOTO 140 - GOTO 90 -C - 120 IF (TYPEBO(J) .EQ. ' LO') GOTO 240 - IF (TYPEBO(J) .EQ. ' MI') GOTO 250 - GOTO 150 -C - 130 IF (TYPEBO(J) .EQ. ' UP') GOTO 260 - IF (TYPEBO(J) .EQ. ' PL') GOTO 200 - GOTO 150 -C - 140 IF (TYPEBO(J) .EQ. ' UP') GOTO 270 - IF (TYPEBO(J) .EQ. ' PL') GOTO 180 -C - 150 WRITE(6,160) TYPEBO(I),TYPEBO(J) - 160 FORMAT('ERROR. Inconsistency between TYPEBO(I) = ', - 1 A3,' and TYPEBO(J) = ',A3) - GOTO 310 -C - 170 CALL FIXVAR(VALBOU(I),COLBOU(I)) - GOTO 290 -C - 180 TYPEBO(J)= BLANK1 - 190 JSTAR=INDCOL(COLBOU(I)) - NFREE= NFREE+1 - VRFREE(NFREE)= ORGIDX(JSTAR) - GOTO 290 -C - 200 TYPEBO(J)= BLANK1 - 210 CALL LOBOVA(VALBOU(I),COLBOU(I)) - GOTO 290 -C - 220 IF (VALBOU(I) .GT. ZEROD) GOTO 230 - CALL CHSIGN(VALBOU(I),COLBOU(I)) - GOTO 210 -C - 230 CALL LOUPVA(VALBOU(I),ZEROD,COLBOU(I)) - GOTO 290 -C - 240 IF (VALBOU(J) .NE. VALBOU(I)) GOTO 245 - TYPEBO(J)= BLANK1 - GOTO 170 - 245 CALL LOUPVA(VALBOU(I),VALBOU(J),COLBOU(I)) - GOTO 280 -C - 250 CALL CHSIGN(VALBOU(I),COLBOU(I)) - GOTO 200 -C - 260 IF (VALBOU(J) .NE. VALBOU(I)) GOTO 265 - TYPEBO(J)= BLANK1 - GOTO 170 - 265 CALL LOUPVA(VALBOU(J),VALBOU(I),COLBOU(I)) - GOTO 280 -C - 270 CALL CHSIGN(VALBOU(J),COLBOU(J)) - CALL LOBOVA(VALBOU(J),COLBOU(J)) -C - 280 TYPEBO(J)= BLANK1 - 290 TYPEBO(I)= BLANK1 -C - 300 CONTINUE -C - 310 RETURN - END -C********************************************************************* -C Function Indcol: associates an index 1, 2, ..n to each column name* -C********************************************************************* - FUNCTION INDCOL(COLX) - INCLUDE (common.h) - CHARACTER*8 COLX -C - DO 10 I=1,NSPCOLS - IF (NAMECL(I) .EQ. COLX) THEN - INDCOL=I - RETURN - ENDIF - 10 CONTINUE -C - WRITE(6,90) COLX - DO 25 I=1,NSPCOLS - WRITE(6,91) NAMECL(I) - 25 CONTINUE - INDCOL= -1 - RETURN - 90 FORMAT(20X,'WARNING. Column index OF ',A8,' was not found.') - 91 FORMAT('NAMECL=',A8) - END -C********************************************************************* -C Function Indrow: associates an index 1, 2, ..n to each row name * -C********************************************************************* - FUNCTION INDROW(ROWX) - INCLUDE (common.h) - CHARACTER*8 ROWX -C - DO 10 I=1,NROWS+1 - IF (NAMERW(I) .EQ. ROWX) THEN - INDROW=I-1 - RETURN - ENDIF - 10 CONTINUE -C - WRITE(6,*) ' ROW:>',ROWX,'',NAMERW(i),'<' - 15 continue - INDROW= -1 - RETURN - END -C********************************************************************* -C Subroutine Sortcl sorts, inside each column, the nonzero coeficients -C in increasing order of row index "JA(I)". * -C********************************************************************* -C - SUBROUTINE SORTCL - INCLUDE (common.h) -C - DOUBLE PRECISION AUX - INTEGER BEGCL1, BEGCL2, ENDCL1, ENDCL2 -C - DO 30 J=1,NSPCOLS - NOZERO= IA(J+1)-IA(J) - IF (NOZERO .LT. 2) GOTO 30 - BEGCL1=IA(J) - ENDCL1=IA(J+1)-2 - ENDCL2=IA(J+1)-1 - DO 20 I1=BEGCL1,ENDCL1 - MIN=I1 - BEGCL2=I1+1 - DO 10 I2=BEGCL2,ENDCL2 - IF (JA(I2) .LT. JA(MIN)) MIN=I2 - 10 CONTINUE - AUX= A(MIN) - IAUX= JA(MIN) - A(MIN)=A(I1) - JA(MIN)=JA(I1) - A(I1)= AUX - JA(I1)=IAUX - 20 CONTINUE - 30 CONTINUE -C - RETURN - END -C -C********************************************************************* -C Subroutine Arwise, stores the coeficient matrix A in a 'rowise' form -C through the arrays JAT, IAT, AT. -C Assumption: within a column in JA, the rows associated with nonzero* -c entries are sorted in crescent order. * -c -C modified to deal with schur complement 6/20/92 -C********************************************************************* -C - SUBROUTINE ARWISE - INCLUDE (common.h) -C - J=1 - DO 30 KROW=1,NROWS - IAT(KROW)= J - DO 20 I1=1, NSPCOLS - DO 10 I2= IA(I1), IA(I1+1)-1 - IF (JA(I2) .EQ. KROW) THEN - JAT(J)= I1 - AT(J)= A(I2) - J= J+1 - GOTO 20 - ELSE - IF (JA(I2) .GT. KROW) GOTO 20 - ENDIF - 10 CONTINUE - 20 CONTINUE - 30 CONTINUE - IAT(NROWS+1)= J -C .............................................Fill in Dense transpose - IF (.NOT. KSCHUR) RETURN - J=1 - DO 90 KROW=1,NROWS - IATD(KROW)= J - DO 80 I1=1, NDNCOLS - DO 70 I2= IAD(I1), IAD(I1+1)-1 - IF (JAD(I2) .EQ. KROW) THEN - JATD(J)= I1 - ATD(J)= AD(I2) - J= J+1 - GOTO 80 - ELSE - IF (JAD(I2) .GT. KROW) GOTO 80 - ENDIF - 70 CONTINUE - 80 CONTINUE - 90 CONTINUE - IATD(NROWS+1)= J -C - RETURN - END ##### sed > DualAffCode/tricks.f <<'#####' 's/^-//' -C ******************************************************************** -C scale: Scales rows, etc. as is described in ob1 manual -C ******************************************************************** - SUBROUTINE SCALE - INCLUDE (common.h) -C - INTEGER*4 IATOAT(MXNZS) -C -C ...... TO operate with rows, I'm going to set up a temporary data -C ...... structure that creates IAT, JAT once, but also creates a -C ...... pointer to the element in A to which the JAT(.)'th element -C ...... corresponds to. -C -C NOTE: The function SCFAC will only work with this data structure. -C .............................................. Form AT. - J=1 - DO 30 KROW=1,NROWS - IAT(KROW)= J - DO 20 I1=1, NSPCOLS - DO 10 I2= IA(I1), IA(I1+1) - 1 - IF (JA(I2) .EQ. KROW) THEN - JAT(J)= I1 - IATOAT(J)=I2 - J= J+1 - GOTO 20 - ELSE - IF (JA(I2) .GT. KROW) GOTO 20 - ENDIF - 10 CONTINUE - 20 CONTINUE - 30 CONTINUE - IAT(NROWS+1)= J -C - DO 50 I=1,NROWS - ROWSF(I)=DBONE - 50 CONTINUE - DO 60 I=1,NSPCOLS - COLSF(I)=DBONE - 60 CONTINUE -C -C ................................... Find each row's GM and divide -C - DO 130 I=1,NROWS - RWGMN1=SCFAC( 2, 1, I, IATOAT) - DO 120 J= IAT(I), IAT(I+1)-1 - A( IATOAT(J) ) = A( IATOAT(J) ) / RWGMN1 - 120 CONTINUE - B(I) = B(I) / RWGMN1 - ROWSF(I)= ROWSF(I) / RWGMN1 - 130 CONTINUE -C ................................... Find each col's GM and divide -C - DO 150 I=1,NSPCOLS - CLGMN1=SCFAC( 1, 1, I, IATOAT) - DO 140 J= IA(I), IA(I+1)-1 - A( J ) = A( J ) / CLGMN1 - 140 CONTINUE - C(I) = C(I) / CLGMN1 - COLSF(I) = COLSF(I) / CLGMN1 - 150 CONTINUE -C ................................... Find each row's GM and divide -C - DO 230 I=1,NROWS - RWGMN2=SCFAC( 2, 1, I, IATOAT) - DO 220 J= IAT(I), IAT(I+1)-1 - A( IATOAT(J) ) = A( IATOAT(J) ) / RWGMN2 - 220 CONTINUE - B(I) = B(I) / RWGMN2 - ROWSF(I)= ROWSF(I) / RWGMN2 - 230 CONTINUE -C ................................... Find each col's GM and divide -C - DO 250 I=1,NSPCOLS - CLGMN2=SCFAC( 1, 1, I, IATOAT) - DO 240 J= IA(I), IA(I+1)-1 - A( J ) = A( J ) / CLGMN2 - 240 CONTINUE - C(I) = C(I) / CLGMN2 - COLSF(I) = COLSF(I) / CLGMN2 - 250 CONTINUE -C ................................... Find each row's LAV and divide -C - DO 330 I=1,NROWS - RWMAXA=SCFAC( 2, 2, I, IATOAT) - DO 320 J= IAT(I), IAT(I+1)-1 - A( IATOAT(J) ) = A( IATOAT(J) ) / RWMAXA - 320 CONTINUE - B(I) = B(I) / RWMAXA - ROWSF(I)= ROWSF(I) / RWMAXA - 330 CONTINUE -C ................................... Find each col's LAV and divide -C - DO 350 I=1,NSPCOLS - CLMAXA=SCFAC( 1, 1, I, IATOAT) - DO 340 J= IA(I), IA(I+1)-1 - A( J ) = A( J ) / CLMAXA - 340 CONTINUE - C(I) = C(I) / CLMAXA - COLSF(I) = COLSF(I) / CLMAXA - 350 CONTINUE -C .................................................... DONE ....... - RETURN - END -C -C **********************************************************8 -C SCFAC(IOPT, JOPT, IRC, IATOAT) : Finds scaling factor of -C IOPT = 1 : Operates on column IRC row or col. -C IOPT = 2 : Operates on row IRC -C JOPT = 1 : Finds Geometric Mean -C JOPT = 2 : Finds largest abs. value. -C Returns requested scaling factor. -C *********************************************************** -C - FUNCTION SCFAC( IOPT, JOPT, IRC, IATOAT) - INCLUDE (common.h) - INTEGER*4 IATOAT(MXNZS) -C - IF(IOPT.EQ.1) THEN - VALMAX = -UBOUND - VALMIN = UBOUND - DO 10 K=IA(IRC), IA(IRC+1)-1 - DTEMP=DABS(A(K)) - IF (DTEMP.LE. EPSILO) GOTO 10 - IF(VALMAX.LT.DTEMP) VALMAX=DTEMP - IF(VALMIN.GT.DTEMP) VALMIN=DTEMP - 10 CONTINUE - ELSE IF(IOPT.EQ.2) THEN - VALMAX = -UBOUND - VALMIN = UBOUND - DO 20 K=IAT(IRC), IAT(IRC+1)-1 - DTEMP=DABS( A(IATOAT(K)) ) - IF (DTEMP.LE. EPSILO) GOTO 20 - IF(VALMAX.LT.DTEMP) VALMAX=DTEMP - IF(VALMIN.GT.DTEMP) VALMIN=DTEMP - 20 CONTINUE - ELSE - WRITE(6,*) ' ERROR: Called SCFAC(',IOPT,JOPT,IRC,')' - ENDIF - IF(JOPT.EQ.1) THEN - SCFAC= DSQRT( VALMIN*VALMAX ) - ELSE IF(JOPT.EQ.2) THEN - SCFAC= VALMAX - ELSE - WRITE(6,*) ' ERROR: Called SCFAC(',IOPT,JOPT,IRC,')' - ENDIF - IF( VALMAX .LT. VALMIN ) SCFAC=DBONE - RETURN - END -C ********************************************************************* -C Subroutine Splitdense: Determines which cols are dense and fills in -C the dense data structures. -C -C Strategy: Separate cols from (IA,JA,A) to a temporary array set -C (IAS, JAS, AS) and permanent array set (IAD, JAD, AD). Then adjust -C NSPCOLS, NDNCOLS, and set IDFLAG -C ********************************************************************* -C - SUBROUTINE SPLITDENSE - INCLUDE (common.h) -C - INTEGER IAS(MXCLS), JAS(MXNZS) - DOUBLE PRECISION AS(MXNZS) -C - NSPC = ZEROI - NDNC = ZEROI - NSPNZ= ZEROI - NDNNZ= ZEROI -C - IF(.NOT.(KSCHUR .OR. KSPLIT)) GOTO 300 - DO 100 ICOL=1, NSPCOLS - IF( (KSCHUR .AND. IA(ICOL+1)-IA(ICOL).GE.IDNTHR) .OR. - 1 (KSPLIT .AND. IDFLAG(ICOL))) THEN -C................................................... Copy into Dense - IDFLAG(ICOL)=.TRUE. - NDNC = NDNC+1 - MPDN2C(NDNC) = ICOL - IAD(NDNC)=NDNNZ+1 - DO 20 J = IA(ICOL), IA(ICOL+1)-1 - NDNNZ=NDNNZ+1 - JAD(NDNNZ)= JA(J) - AD(NDNNZ) = A(J) - 20 CONTINUE - ELSE -C................................................... Copy into SPARSE - if(idflag(icol)) then - write(0,*) ' AAAAH: idflag s/b false.' - stop - endif - IDFLAG(ICOL)=.FALSE. - NSPC = NSPC+1 - MPSP2C(NSPC) = ICOL - IAS(NSPC)=NSPNZ+1 - DO 30 J = IA(ICOL), IA(ICOL+1)-1 - NSPNZ=NSPNZ+1 - JAS(NSPNZ)= JA(J) - AS(NSPNZ) = A(J) - 30 CONTINUE - ENDIF - 100 CONTINUE - IAS(NSPC+1)=NSPNZ+1 -C................................................... Copy AS back into A. - DO 220 ICOL=1, NSPC - IA(ICOL)=IAS(ICOL) - DO 210 JROW= IAS(ICOL), IAS(ICOL+1)-1 - JA(JROW)= JAS(JROW) - A(JROW)= AS(JROW) - 210 CONTINUE - 220 CONTINUE - IAD(NDNC+1)=NDNNZ+1 - IA(NSPC+1) =NSPNZ+1 - NSPCOLS =NSPC -C................................................. Come here if no schur comp. - 300 NDNCOLS =NDNC - NCOLS =NSPCOLS+NDNCOLS - 777 format(3(2x,i4,2x,A8,3x,I4)) -C - RETURN - END -C **************************************************************** -C Subroutine ColCount : counts and prints the no. of nonzeroes in -C each column. -C By dfholmes, 6/20 -C ***************************************************************** - SUBROUTINE COLCOUNT - INCLUDE (common.h) -C - WRITE(6,900) - DO 10 I=1,NSPCOLS - ITEMP=ZEROI - IRWCNT=IA(I+1)-IA(I) - DO 5 J=IA(I),IA(I+1)-2 - ITEMP=ITEMP+JA(J+1)-JA(J)+1 - 5 CONTINUE - IDFLAG(I)=.FALSE. - DTEMP = DFLOAT(ITEMP)/DFLOAT(IRWCNT) - DTEMP2 = DFLOAT(ITEMP)/DFLOAT(NROWS) - IF(KSPLIT) THEN - IF( DTEMP .GE. DFLOAT(IDNTHR) ) THEN - IDFLAG(I)=.TRUE. - ENDIF - ELSE - WRITE(6,910) I, NAMECL(I), IRWCNT, DTEMP, DTEMP2 - ENDIF - 10 CONTINUE -C - 900 FORMAT('---------- Column Nonzero Count ---------------'/ - 1 ' Avg Rows Col. '/ - 1 ' Column Name Nonzeroes bet. nonz Density'/ - 2 ' -----------------------------------------------') - 910 FORMAT(1X,I5, 2X, A8, 2X,I8,2X,F11.6,2X,F8.6) - RETURN - END -C **************************************************************** -C Subroutine FINDDG : calculates the diagonal entries in D for AD^2 A^T. -C Modify this routine to try the primal-method or the primal-dual method. -C By dfholmes, 6/20 -C ***************************************************************** - SUBROUTINE FINDDG - INCLUDE (common.h) -C - IF (IALG .EQ. IDUALA) THEN - DO 10 I=1,NSPCOLS+NDNCOLS - D(I) = DBONE/(SLACK(I)*SLACK(I)) - DROOT(I) = DBONE/SLACK(I) - 10 CONTINUE - ENDIF -C - RETURN - END ##### sed > DualAffCode/wrapup.f <<'#####' 's/^-//' -C ************************************************************************ -C Subroutine Dualiz, selects which problem to solve, the dual C or primal, -C by choosing the one which has the least number of constraints. If the -C algorithm will be applied to the primal then it exchanges A by At using -C the arrays IA, JA and A. It also exchanges -C and B. -C ************************************************************************ - SUBROUTINE DUALIZ - INCLUDE (common.h) - INTEGER IAUX(MXNZS) - DOUBLE PRECISION AUX(MXNZS) -C - IF (PROBLM .EQ. 'DUAL') GOTO 120 - IF (PROBLM .EQ. 'PRML') GOTO 10 - IF (NROWS .LE. NCOLS) THEN - PROBLM= 'DUAL' - GOTO 120 - ELSE - PROBLM= 'PRML' - END IF -C - 10 IEND1= NCOLS+1 - DO 20 I=1,IEND1 - IAUX(I)= IA(I) - 20 CONTINUE -C - IEND2= NROWS+1 - DO 30 I= 1,IEND2 - IA(I)= IAT(I) - 30 CONTINUE -C - DO 40 I= 1, IEND1 - IAT(I)= IAUX(I) - 40 CONTINUE -C - IEND1= IAT(NCOLS+1)-1 - DO 50 I= 1, IEND1 - IAUX(I)= JA(I) - AUX(I)= A(I) - 50 CONTINUE -C - IEND2=IA(NROWS+1)-1 - DO 60 I=1,IEND2 - JA(I)= JAT(I) - A(I)= AT(I) - 60 CONTINUE -C - DO 70 I=1,IEND1 - JAT(I)= IAUX(I) - AT(I)= AUX(I) - 70 CONTINUE -C - DO 80 I=1,NROWS - AUX(I)= B(I) - 80 CONTINUE - DO 90 I=1,NCOLS - B(I)= -C(I) - 90 CONTINUE - IEND1= NCOLS+1 - DO 100 I=IEND1,NROWS - B(I)= ZEROD - 100 CONTINUE - DO 110 I=1,NROWS - C(I)=AUX(I) - 110 CONTINUE -C - IAUX(1)= NROWS - NROWS= NCOLS - NCOLS= IAUX(1) -C - 120 RETURN - END -C -C ************************************************************* -C * * -C * SUBROUTINE PRINT : used to print out the final results, * -C * * -C ************************************************************* -C - SUBROUTINE PRINT - INCLUDE (common.h) -C - DOUBLE PRECISION APRINT(20,20) -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 - IF (FEASBL .EQ. 'STUCK-3 ') THEN - WRITE(6,860) 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.0 -C 120 CONTINUE -C 130 CONTINUE -c ---------------------------- must change this code for schur compl. -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,609) (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 - ITIM(7)=MCLOCK() - WRITE(6,760) - WRITE(6,730) ZLOWER - WRITE(6,740) ZUPPER - IF (ARTFCL .LE. 0.0) 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) NDNCOLS - WRITE(6,491) NSPCOLS - WRITE(6,492) NCOLS - WRITE(6,493) NZEROP - WRITE(6,500) ZCONST - WRITE(6,530) NRWDRP - WRITE(6,540) NFIX - WRITE(6,570) NSTFRE - WRITE(6,550) PROBLM - WRITE(6,560) PHASE - WRITE(6,580) DFLOAT(ITIM(3)-ITIM(2))/100.0D0 - WRITE(6,590) DFLOAT(ITIM(2)-ITIM(1))/100.0D0 - WRITE(6,600) DFLOAT(ITIM(4)-ITIM(3))/100.0D0 - WRITE(6,601) DFLOAT(ITIM(10))/100.0D0 - WRITE(6,602) DFLOAT(ITIM(11))/100.0D0 - WRITE(6,603) DFLOAT(ITIM(12))/100.0D0 - WRITE(6,610) - WRITE(6,620) DFLOAT(ITIM(5)-ITIM(4))/100.0D0, - 1 DFLOAT(ITIM(6)-ITIM(5))/100.0D0 - WRITE(6,750) ITERA1, ITERA2 - WRITE(6,650) DFLOAT(ITIM(7)-ITIM(6))/100.0D0 - WRITE(6,630) DFLOAT(ITIM(7)-ITIM(1))/100.0D0 -C - 480 FORMAT(10X,'# of ORIGINAL constraints ..',I14) - 490 FORMAT(10X,'# of DENSE columns ......',I14) - 491 FORMAT(10X,'# of SPARSE columns ......',I14) - 492 FORMAT(10X,'# of ORIGINAL columns ......',I14) - 493 FORMAT(10X,'# Zero Pivot Rows ..........',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(10X,'Search Direction time ......',F14.5) - 602 FORMAT(10X,'Step Size time .............',F14.5) - 603 FORMAT(10X,'Update time ................',F14.5) - 609 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= ',E30.20) - 740 FORMAT(10X,'Upper Bound on Z ...... ZUPPER= ',E30.20/) - 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:',//) - 860 FORMAT(/,10X,'WARNING: Y(ARTF) BQ_Decomp/blkio.f <<'#####' 's/^-//' -C -C*************************************************************** -C * -C BLKIO keeps track (stores and retrieves) of factorizations * -C of S. Called by FINDDY * -C * -C*************************************************************** -C - SUBROUTINE BLKIO ( IOPT, iblk, S ) - INCLUDE (sizes.h) -C -C*************************************************************** -C -C IBLK = 1...ISCEN -C -C IOPT=1 : WRITE NOTE: Program requires blocks to be saved -C IOPT=2 : READ in increasing order (not nec. for par. imp.) -C -C ISAVE( 1 = NVARS RSAVE() = SVALUE -C 2 = STAGE -C 3 = MXUSED IDSAVE(IBLK) -> s(1) FOR IBLK -C 4-53 = ICPAD IDSAVE(IBLK+1)-1 -> S(LST) FOR IBLK -C 54-101= MAPPTR -C -C*************************************************************** -C - INTEGER I , ICPAD , IERRA , IPRNTE, IPRNTS, - 1 IUNIT , MAPPTR, MAXINT, MAXSA , MSGLVA, - 1 MXUSED, NVARS , STAGE, ISAVE, IDSAVE, IOPT, - 1 LEN - REAL MCHEPS, RATIOL, RATIOS, SVALUE, CLOCK, RSAVE - DOUBLE PRECISION S(SIZEMAT), DSAVE -C -C*************************************************************** -C - COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, - 1 MCHEPS, CLOCK - COMMON /SPAUSR/ MSGLVA, IERRA , MAXSA , NVARS - 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*************************************************************** -C - IF (IOPT.EQ.1) THEN -C - ISAVE(IBLK,1)=NVARS - ISAVE(IBLK,2)=STAGE - ISAVE(IBLK,3)=MXUSED -C - DO 10 I1=1,50 - RSAVE(IBLK,I1) = SVALUE(I1) - ISAVE(IBLK,3+I1) = ICPAD(I1) - IF(I1.LE.48) ISAVE(IBLK,53+I1)=MAPPTR(I1) - 10 CONTINUE -C - IF (MXUSED.NE.0) THEN - DO 20 I=1,MXUSED - DSAVE(I+IDSAVE(IBLK)-1)=S(I) - 20 CONTINUE - IDSAVE(IBLK+1)=IDSAVE(IBLK)+MXUSED - ENDIF -C - ELSE IF (IOPT.EQ.2) THEN -C - NVARS =ISAVE(IBLK,1) - STAGE =ISAVE(IBLK,2) - MXUSED=ISAVE(IBLK,3) -C - DO 30 I1=1,50 - SVALUE(I1) = RSAVE(IBLK,I1) - ICPAD(I1) = ISAVE(IBLK,3+I1) - IF(I1.LE.48) MAPPTR(I1) = ISAVE(IBLK,53+I1) - 30 CONTINUE -C - LEN = IFIX(FLOAT(MAXSA)*RATIOL + 0.01) - CALL ZEROIL ( LEN, S(1) ) -C - IF (MXUSED.NE.0 ) THEN - DO 40 I=1,MXUSED - S(I)=DSAVE(I+IDSAVE(IBLK)-1) - 40 CONTINUE - ENDIF -C - ELSE -C - WRITE(6,*) ' .... Error, BLKIO: iopt not in {1,2} ......' - STOP - ENDIF -C - RETURN - END ##### sed > BQ_Decomp/finddy.f <<'#####' 's/^-//' -C ***************************************************************** -C * * -C * SUB: FINDDY: THIS is the meat of the matter. * -C * This subroutine uses Birge and Qi's ddevelopment * -C * to solve the system AD2ATdy=b. * -C * The output of the subroutine is the direction of * -C * improvement dy. * -C * * -C * Called by Main Program, subroutines called: (Sparspak) * -C ***************************************************************** -C - SUBROUTINE FINDDY(ITERTN,U,DY) - INCLUDE (sizes.h) -C - DOUBLE PRECISION G1(SIZET,SIZET),ytmp1(SIZEVEC),ytmp2(SIZEVEC) - DOUBLE PRECISION phat1(SIZEVEC), ytmp0(SIZEVEC),qq(SIZEVEC) - DOUBLE PRECISION DY(SIZEVEC), A,B,C,ZCONST,X,Y,ytmp3(SIZEVEC) - DOUBLE PRECISION AT,AAT,DIAG,ZEROD,SLACK,DTEMP,G2(SIZET,SIZET) - DOUBLE PRECISION S(SIZEMAT), DSUM, U(SIZEVEC), DSAVE -C - INTEGER IA,JA,IAT,JAT,IAAT,JAAT,IDIAG,NCOLS,NROWS,ZEROI,MXITNS - INTEGER ITERA1,ITERA2,INDX(SIZEBLK),INDX2(SIZEBLK),ISAVE,IDSAVE - INTEGER ITERTN -C - REAL RSAVE -C - CHARACTER*3 BLANK1 - CHARACTER*8 BLANK2 -C - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, - 1 STUPTM, PHS1TM, PHS2TM, MXITNS - 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), - 2 JAAT(SIZEAAT) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /SPAUSR/ MSGLVA,IERRA,MAXSA, NEQNS - 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 /SPSAVE/ ISAVE(MAXBLK,103),RSAVE(MAXBLK,50),IDSAVE(MAXBLK), - 1 DSAVE(SIZEDSAVE) -C -C .................. Initialize ..................... -C - MAXSA=SIZEMAT - IDSAVE(1)=1 - MSGLVA=1 -C -C .................. SOLVE Sp=b ..................... -C -C Solve P(0) -C -C------------------------------------------------------------------------- -c write(6,*) ' ============== F I N D D Y =======================' -C------------------------------------------------------------------------- -c if(ITERTN.eq.3) MSGLVA=4 - DO 300 I=1,IBLSZ(1,2) - DY(I)=U(I) - 300 CONTINUE - IRWOFS=IBLSZ(1,2) -C /* Form A0tAo */ - DO 314 I1=IBLA(1),IBLA(2)-1 - DO 313 I2=IBLA(1),IBLA(2)-1 - DSUM=ZEROD - DO 312 I3=IA(I1),IA(I1+1)-1 - DO 311 I4=IA(I2),IA(I2+1)-1 - IF(JA(I3).eq.JA(I4)) DSUM=DSUM+A(I3)*A(I4) - 311 CONTINUE - 312 CONTINUE - G1(I1,I2)=DSUM - 313 CONTINUE - 314 CONTINUE -C -C -C Solve P(I) (directly into dy) -C -C Note, to save computation, I'm also going to compute G1 siumltaneously -C -C /* Input Structure */ - DO 450 IBLK=1,ISCEN - I0=1+ISCEN+IBLK -C !!!! Check to make sure this zeroes A,b !!! - IF(ITERTN.eq.1) THEN - CALL IJBEGN - DO 330 I1=IBLAAT(I0),IBLAAT(I0+1)-1 - ICOL=I1-IBLAAT(I0)+1 - CALL INIJ(ICOL,ICOL,S) - DO 325 I2=IAAT(I1),IAAT(I1+1)-1 - CALL INIJ(JAAT(I2),ICOL,S) - 325 CONTINUE - 330 CONTINUE - CALL IJEND(S) -c MSGLVA=1 - CALL ORDRA1(S) SPK01 -c MSGLVA=0 - ELSE - CALL BLKIO(2,IBLK,S) - ENDIF -C /* Input Numerics */ -C /* Assumes Symmetric */ - DO 340 I1=IBLAAT(I0),IBLAAT(I0+1)-1 - DO 335 I2=IAAT(I1),IAAT(I1+1)-1 - ICOL=I1-IBLAAT(I0)+1 - JTEMP=JAAT(I2) -C IF(JTEMP.LT.I1) GOTO 335 - CALL INAIJ1(JTEMP,ICOL,AAT(I2),S) SPK02 - 335 CONTINUE - 340 CONTINUE - DO 350 I1=IDIAG(I0),IDIAG(I0+1)-1 - I2=I1-IDIAG(I0)+1 - CALL INAIJ1(I2,I2,DIAG(I1),S) SPK03 - 350 CONTINUE -C /* Solve for p*/ - DO 360 I1=1,IBLSZ(I0,2) - CALL INBI(I1,U(I1+IRWOFS),S) - 360 CONTINUE - CALL SOLVE1(S) SPK04 - DO 370 I1=1,IBLSZ(I0,2) - dy(I1+IRWOFS)=S(I1) - 370 CONTINUE -c write(8,949) (i+irwofs,dy(irwofs+i),i=1,iblsz(i0,2)) -c949 format(2(' initdy(',i3,')=',E13.7)) -c - IRWOFS=IRWOFS+IBLSZ(I0,2) -C /* Load Col A(i) */ -C /* YTMP1=Inv(S)*A */ -C /* YTMP2=At*YTMP1 */ - ngtotl1=0 - DO 400 I1=IBLA(IBLK+1),IBLA(IBLK+2)-1 - ICOLI=(I1-IBLA(IBLK+1)+1)+IACLOF - DO 380 I2=IA(I1),IA(I1+1)-1 - CALL INBI(JA(I2),A(I2),S) - 380 CONTINUE - CALL SOLVE1(S) - DO 385 I2=1,IBLSZ(IBLK+1,2) - YTMP1(I2)=S(I2) - 385 CONTINUE -C /* Can't call MULATY, the cols are offset */ -C - DO 387 I2=IBLA(IBLK+1),IBLA(IBLK+2)-1 - DSUM=ZEROD - DO 386 I3=IA(I2),IA(I2+1)-1 - DSUM=DSUM+A(I3)*YTMP1(JA(I3)) - 386 CONTINUE - JROWJ=(I2-IBLA(IBLK+1)+1)+IACLOF - G1(JROWJ,ICOLI)=G1(JROWJ,ICOLI)+DSUM - ngtotl1=ngtotl1+1 - 387 CONTINUE - 400 CONTINUE -C /* Save Ch.Fac of S(i) */ - CALL BLKIO(1,IBLK,S) -C /* End Main Loop */ - 450 CONTINUE -C /* Complete G1 - DO 480 I1=1,IBLSZ(1,1) - G1(I1,I1)=G1(I1,I1)+(SLACK(I1)*SLACK(I1)) - 480 CONTINUE -c -C------------------------------------------------------------------------- -c write(8,*) ' IT: ',ITERTN,'Final G1, ngtotls=',ngtotl1 -c do 482 i=1,iblsz(1,1) -c do 481 j=1,iblsz(1,1) -c if(dabs(g1(i,j)).GE.1.0D-03) write(8,950) itertn, i,j,g1(i,j) -c481 continue -c482 continue - 950 format(2(' IT:',I2,'g1(',i3,',',i3,')=',E13.7)) - 951 format(2(' IT:',I2,'g2(',i3,',',i3,')=',E13.7)) - 952 format(3(' qq(',I4,')=',E13.7)) - 953 format(3(' u(',I4,')=',E13.7)) - 954 format(3(' y(',I4,')=',E13.7)) - 955 format(3(' t(',I4,')=',E13.7)) - 956 format(3(' ph(',I4,')=',E13.7)) -C------------------------------------------------------------------------- - -C -C .................. SOLVE Gq=Vtp ................... -C -C /* clr phat1 \in \Re^NACOLS */ -C DO 500 I1=1,IACLOF+IBLSZ(1,1) -C PHAT1(I1)=ZEROD -C500 CONTINUE -C -C /* A0tp(0) */ - DO 508 I1=IBLA(1),IBLA(2)-1 - DSUM=ZEROD - DO 507 I2=IA(I1),IA(I1+1)-1 - DSUM=DSUM+A(I2)*DY(JA(I2)) - 507 CONTINUE - PHAT1(I1-IBLA(1)+1)=DSUM - 508 CONTINUE -C /* IROWOFS=NAROWS */ - IRWOFS=IBLSZ(1,2) -C /* Add other Aitp(i) */ - DO 520 IBLK=1,ISCEN - DO 513 I1=IBLA(IBLK+1),IBLA(IBLK+2)-1 - DSUM=ZEROD - DO 510 I2=IA(I1),IA(I1+1)-1 - DSUM=DSUM+A(I2)*DY(IRWOFS+JA(I2)) - 510 CONTINUE - I3=I1-IBLA(IBLK+1)+1+IACLOF - PHAT1(I3)=PHAT1(I3)+DSUM -C - 513 CONTINUE - IRWOFS=IRWOFS+IBLSZ(IBLK+1,2) - 520 CONTINUE -C /* Solve G1x=phat1 */ - DO 547 I1=1,IBLSZ(1,1) - YTMP2(I1)=PHAT1(I1) - 547 CONTINUE -C - CALL LUDCMP(G1,IBLSZ(1,1),SIZET,INDX,DTEMP) -C - CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP2) -C /* Form ytmp0=-Ax */ - DO 555 I1=IBLAT(1),IBLAT(2)-1 - DSUM=ZEROD - DO 550 I2=IAT(I1),IAT(I1+1)-1 - DSUM=DSUM+AT(I2)*YTMP2(JAT(I2)) - 550 CONTINUE - I3=I1-IBLAT(1)+1 - IF(I3.GT.IBLSZ(1,2)) WRITE(6,*) 'ERROR2 IN GQ=UTP' - YTMP0(I3)=+DY(I3)-DSUM - 555 CONTINUE -C /* Form G2, G1 already factored*/ - DO 580 I1=IBLAT(1),IBLAT(2)-1 -C /* S/B NAROWS */ - DO 570 I9=1,iblsz(1,1) - YTMP1(I9)=0.D0 - 570 CONTINUE -C /* solve g1x=ytmp1 */ - DO 571 I2=IAT(I1),IAT(I1+1)-1 - YTMP1(JAT(I2))=AT(I2) - 571 CONTINUE -C - CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP1) -C -C /* -A* */ - DO 575 I2=IBLAT(1),IBLAT(2)-1 - I2A=I2-IBLAT(1)+1 - DSUM=ZEROD - DO 574 I3=IAT(I2),IAT(I2+1)-1 - DSUM=DSUM-AT(I3)*YTMP1(JAT(I3)) - 574 CONTINUE - G2(I2A,I1)=DSUM - 575 CONTINUE - 580 CONTINUE -C /* Now solve g2q2=y for q2 */ -C /* First, input G2 into SPA */ -C - CALL LUDCMP(G2,IBLSZ(1,2),SIZET,INDX2,DTEMP) -C - CALL LUBKSB(G2,IBLSZ(1,2),SIZET,INDX2,YTMP0) -C------------------------------------------------------------------------- -c write(8,*) ' IT: ',ITERTN,'G2.............' -c do 582 i=1,iblsz(1,2) -c do 581 j=1,iblsz(1,2) -c write(8,951) itertn, i,j,g2(i,j) -c581 continue -c582 continue -c------------------------------------------------------------------------- */ - -C - DO 605 I1=1,IBLSZ(1,2) - QQ(I1+IBLSZ(1,1))=YTMP0(I1) - 605 CONTINUE -C -C /* Form phat1-At.q2 */ -C /* First, clear ytmp2 */ -C DO 607 I1=1,3000 -C YTMP3(I1)=ZEROD -C607 CONTINUE -C - DO 615 I1=IBLA(1),IBLA(2)-1 - I3=I1-IBLA(1)+1 - YTMP3(I3)=PHAT1(I3) - DO 610 I2=IA(I1),IA(I1+1)-1 - YTMP3(I3)=YTMP3(I3)-A(I2)*QQ(IBLSZ(1,1)+JA(I2)) - 610 CONTINUE - 615 CONTINUE -C /* SOLVE G1q=ytmp3(.) */ - CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP3) -C - DO 626 I1=1,IBLSZ(1,1) - QQ(I1)=YTMP3(I1) - 626 CONTINUE -c------------------------------------------------------------------------- -C /* NOW qq is complete */ -C .................. SOLVE Sr=Uq .................... -C /* Let dy=dy-r, form r(0) */ - DO 635 I1=IBLAT(1),IBLAT(2)-1 - DSUM=ZEROD - DO 630 I2=IAT(I1),IAT(I1+1)-1 - DSUM=DSUM+AT(I2)*QQ(JAT(I2)) - 630 CONTINUE - I3=I1-IBLAT(1)+1 - IF(I3.GT.IBLSZ(1,2)) WRITE(6,*) 'ERROR2 IN Sr=Uq' - DY(I3)=DY(I3)-QQ(I3+IBLSZ(1,1))-DSUM - 635 CONTINUE - IRWOFS=IBLSZ(1,2) -C /* Form r(i), or equiv */ - DO 660 I1=1,ISCEN -C /* Reload Ch.Fac. of S */ - CALL BLKIO(2,I1,S) - DO 645 I2=IBLAT(I1+1),IBLAT(I1+2)-1 - I2A=I2-IBLAT(I1+1)+1 - DSUM=ZEROD - DO 640 I3=IAT(I2),IAT(I2+1)-1 - DSUM=DSUM+AT(I3)*QQ(JAT(I3)+IACLOF) - 640 CONTINUE - CALL INBI(I2A,DSUM,S) - 645 CONTINUE - CALL SOLVE1(S) SPK06 - DO 650 I2=IRWOFS+1,IRWOFS+IBLSZ(IBLK+1,2) - I3=I2-IRWOFS - DY(I2)=DY(I2)-S(I3) - 650 CONTINUE - IRWOFS=IRWOFS+IBLSZ(I1+1,2) - 660 CONTINUE -C------------------------------------------------------------------------- -c write(8,*) 'ITERTN=',ITERTN -c write(8,903) (i,dy(i),i=1,irwofs+1) -c903 format(2(' dyp(',i3,')=',E13.7)) -c904 format(8(' qq(',i3,')=',F9.5)) -C------------------------------------------------------------------------- - -C /* DONE with dy (EX. PH1) */ - RETURN - END ##### sed > BQ_Decomp/finddylin.f <<'#####' 's/^-//' -C ***************************************************************** -C * * -C * SUB: FINDDY: THIS is the meat of the matter. * -C * This subroutine uses Birge and Qi's ddevelopment * -C * to solve the system AD2ATdy=b. * -C * The output of the subroutine is the direction of * -C * improvement dy. * -C * * -C * Called by Main Program, subroutines called: (Sparspak) * -C ***************************************************************** -C - SUBROUTINE FINDDY(ITERTN,U,DY,RCOND) - INCLUDE (sizes.h) -C - DOUBLE PRECISION G1(SIZET,SIZET),ytmp1(SIZEVEC),ytmp2(SIZEVEC) - DOUBLE PRECISION phat1(SIZEVEC), ytmp0(SIZEVEC),qq(SIZEVEC) - DOUBLE PRECISION DY(SIZEVEC), A,B,C,ZCONST,X,Y,ytmp3(SIZEVEC) - DOUBLE PRECISION AT,AAT,DIAG,ZEROD,SLACK,DTEMP,G2(SIZET,SIZET) - DOUBLE PRECISION S(SIZEMAT), DSUM, U(SIZEVEC), DSAVE - DOUBLE PRECISION Z(SIZET), RCOND(4) -C - INTEGER IA,JA,IAT,JAT,IAAT,JAAT,IDIAG,NCOLS,NROWS,ZEROI,MXITNS - INTEGER ITERA1,ITERA2,INDX(SIZEBLK),INDX2(SIZEBLK),ISAVE,IDSAVE - INTEGER ITERTN, KPVT(SIZET) -C - REAL RSAVE -C - CHARACTER*3 BLANK1 - CHARACTER*8 BLANK2 -C - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /COUNT/ ITERA1, ITERA2, PREPTM, READTM, - 1 STUPTM, PHS1TM, PHS2TM, MXITNS - 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), - 2 JAAT(SIZEAAT) - COMMON /BODY/ SLACK(SIZEVEC) - COMMON /BLANK/ ZEROD, ZEROI, BLANK1, BLANK2 - COMMON /SPAUSR/ MSGLVA,IERRA,MAXSA, NEQNS - 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 /SPSAVE/ ISAVE(MAXBLK,103),RSAVE(MAXBLK,50),IDSAVE(MAXBLK), - 1 DSAVE(SIZEDSAVE) -C -C .................. Initialize ..................... -C - MAXSA=SIZEMAT - IDSAVE(1)=1 - MSGLVA=1 -C -C .................. Solve Sp=b ..................... -C -C Solve P(0) -C -C------------------------------------------------------------------------- -c write(6,*) ' ============== F I N D D Y =======================' -C------------------------------------------------------------------------- -c if(ITERTN.eq.3) MSGLVA=4 - DO 300 I=1,IBLSZ(1,2) - DY(I)=U(I) - 300 CONTINUE - IRWOFS=IBLSZ(1,2) -C /* Form A0tAo */ - DO 314 I1=IBLA(1),IBLA(2)-1 - DO 313 I2=IBLA(1),IBLA(2)-1 - DSUM=ZEROD - DO 312 I3=IA(I1),IA(I1+1)-1 - DO 311 I4=IA(I2),IA(I2+1)-1 - IF(JA(I3).eq.JA(I4)) DSUM=DSUM+A(I3)*A(I4) - 311 CONTINUE - 312 CONTINUE - G1(I1,I2)=DSUM - 313 CONTINUE - 314 CONTINUE -C -C -C Solve P(I) (directly into dy) -C -C Note, to save computation, I'm also going to compute G1 siumltaneously -C -C /* Input Structure */ - DO 450 IBLK=1,ISCEN - I0=1+ISCEN+IBLK -C !!!! Check to make sure this zeroes A,b !!! - IF(ITERTN.eq.1) THEN - CALL IJBEGN - DO 330 I1=IBLAAT(I0),IBLAAT(I0+1)-1 - ICOL=I1-IBLAAT(I0)+1 - CALL INIJ(ICOL,ICOL,S) - DO 325 I2=IAAT(I1),IAAT(I1+1)-1 - CALL INIJ(JAAT(I2),ICOL,S) - 325 CONTINUE - 330 CONTINUE - CALL IJEND(S) -c MSGLVA=1 - CALL ORDRA5(S) SPK01 -c MSGLVA=0 - ELSE - CALL BLKIO(2,IBLK,S) - ENDIF -C /* Input Numerics */ -C /* Assumes Symmetric */ - DO 340 I1=IBLAAT(I0),IBLAAT(I0+1)-1 - DO 335 I2=IAAT(I1),IAAT(I1+1)-1 - ICOL=I1-IBLAAT(I0)+1 - JTEMP=JAAT(I2) -C IF(JTEMP.LT.I1) GOTO 335 - CALL INAIJ5(JTEMP,ICOL,AAT(I2),S) SPK02 - 335 CONTINUE - 340 CONTINUE - DO 350 I1=IDIAG(I0),IDIAG(I0+1)-1 - I2=I1-IDIAG(I0)+1 - CALL INAIJ5(I2,I2,DIAG(I1),S) SPK03 - 350 CONTINUE -C /* Solve for p*/ - DO 360 I1=1,IBLSZ(I0,2) - CALL INBI(I1,U(I1+IRWOFS),S) - 360 CONTINUE - CALL SOLVE5(S) SPK04 - RCOND(3)=-1.0D0 - MSGLVA=4 -c CALL EREST5(RCOND(3),S) SPK042 - MSGLVA=1 -c write(6,*) 'CEst: it: ',itertn,' blk: ',iblk,' rcond: ',rcond(3) - DO 370 I1=1,IBLSZ(I0,2) - dy(I1+IRWOFS)=S(I1) - 370 CONTINUE -c write(8,949) (i+irwofs,dy(irwofs+i),i=1,iblsz(i0,2)) -c949 format(2(' initdy(',i3,')=',E13.7)) -c - IRWOFS=IRWOFS+IBLSZ(I0,2) -C /* Load Col A(i) */ -C /* YTMP1=Inv(S)*A */ -C /* YTMP2=At*YTMP1 */ - ngtotl1=0 - DO 400 I1=IBLA(IBLK+1),IBLA(IBLK+2)-1 - ICOLI=(I1-IBLA(IBLK+1)+1)+IACLOF - DO 380 I2=IA(I1),IA(I1+1)-1 - CALL INBI(JA(I2),A(I2),S) - 380 CONTINUE - CALL SOLVE5(S) - DO 385 I2=1,IBLSZ(IBLK+1,2) - YTMP1(I2)=S(I2) - 385 CONTINUE -C /* Can't call MULATY, the cols are offset */ -C - DO 387 I2=IBLA(IBLK+1),IBLA(IBLK+2)-1 - DSUM=ZEROD - DO 386 I3=IA(I2),IA(I2+1)-1 - DSUM=DSUM+A(I3)*YTMP1(JA(I3)) - 386 CONTINUE - JROWJ=(I2-IBLA(IBLK+1)+1)+IACLOF - G1(JROWJ,ICOLI)=G1(JROWJ,ICOLI)+DSUM - ngtotl1=ngtotl1+1 - 387 CONTINUE - 400 CONTINUE -C /* Save Ch.Fac of S(i) */ - CALL BLKIO(1,IBLK,S) -C /* End Main Loop */ - 450 CONTINUE -C /* Complete G1 - DO 480 I1=1,IBLSZ(1,1) - G1(I1,I1)=G1(I1,I1)+(SLACK(I1)*SLACK(I1)) - 480 CONTINUE -c -C------------------------------------------------------------------------- -c write(8,*) ' IT: ',ITERTN,'Final G1, ngtotls=',ngtotl1 -c do 482 i=1,iblsz(1,1) -c do 481 j=1,iblsz(1,1) -c if(dabs(g1(i,j)).GE.1.0D-03) write(8,950) itertn, i,j,g1(i,j) -c481 continue -c482 continue - 950 format(2(' IT:',I2,'g1(',i3,',',i3,')=',E13.7)) - 951 format(2(' IT:',I2,'g2(',i3,',',i3,')=',E13.7)) - 952 format(3(' qq(',I4,')=',E13.7)) - 953 format(3(' u(',I4,')=',E13.7)) - 954 format(3(' y(',I4,')=',E13.7)) - 955 format(3(' t(',I4,')=',E13.7)) - 956 format(3(' ph(',I4,')=',E13.7)) -C------------------------------------------------------------------------- - -C -C .................. SOLVE Gq=Vtp ................... -C -C /* clr phat1 \in \Re^NACOLS */ -C DO 500 I1=1,IACLOF+IBLSZ(1,1) -C PHAT1(I1)=ZEROD -C500 CONTINUE -C -C /* A0tp(0) */ - DO 508 I1=IBLA(1),IBLA(2)-1 - DSUM=ZEROD - DO 507 I2=IA(I1),IA(I1+1)-1 - DSUM=DSUM+A(I2)*DY(JA(I2)) - 507 CONTINUE - PHAT1(I1-IBLA(1)+1)=DSUM - 508 CONTINUE -C /* IROWOFS=NAROWS */ - IRWOFS=IBLSZ(1,2) -C /* Add other Aitp(i) */ - DO 520 IBLK=1,ISCEN - DO 513 I1=IBLA(IBLK+1),IBLA(IBLK+2)-1 - DSUM=ZEROD - DO 510 I2=IA(I1),IA(I1+1)-1 - DSUM=DSUM+A(I2)*DY(IRWOFS+JA(I2)) - 510 CONTINUE - I3=I1-IBLA(IBLK+1)+1+IACLOF - PHAT1(I3)=PHAT1(I3)+DSUM -C - 513 CONTINUE - IRWOFS=IRWOFS+IBLSZ(IBLK+1,2) - 520 CONTINUE -C /* Solve G1x=phat1 */ - DO 547 I1=1,IBLSZ(1,1) - YTMP2(I1)=PHAT1(I1) - 547 CONTINUE -C - RCOND(1)=-1.0D0 - CALL DPOCO(G1,SIZET,IBLSZ(1,1),RCOND(1),Z,INFO) - IF(INFO.NE.0) WRITE(0,*) ' AFTER G1 FACTORIZATION: INFO=',INFO -c write(0,*) ' ITERTN: ',ITERTN,' G1 Cond: ',RCOND(1) - CALL DPOSL(G1,SIZET,IBLSZ(1,1),YTMP2) -C -C CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP2) -C /* Form ytmp0=-Ax */ - DO 555 I1=IBLAT(1),IBLAT(2)-1 - DSUM=ZEROD - DO 550 I2=IAT(I1),IAT(I1+1)-1 - DSUM=DSUM+AT(I2)*YTMP2(JAT(I2)) - 550 CONTINUE - I3=I1-IBLAT(1)+1 - IF(I3.GT.IBLSZ(1,2)) WRITE(6,*) 'ERROR2 IN GQ=UTP' - YTMP0(I3)=+DY(I3)-DSUM - 555 CONTINUE -C /* Form G2, G1 already factored*/ - DO 580 I1=IBLAT(1),IBLAT(2)-1 -C /* S/B NAROWS */ - DO 570 I9=1,iblsz(1,1) - YTMP1(I9)=0.D0 - 570 CONTINUE -C /* solve g1x=ytmp1 */ - DO 571 I2=IAT(I1),IAT(I1+1)-1 - YTMP1(JAT(I2))=AT(I2) - 571 CONTINUE -C - CALL DPOSL(G1,SIZET,IBLSZ(1,1),YTMP1) -C CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP1) -C -C /* -A* */ - DO 575 I2=IBLAT(1),IBLAT(2)-1 - I2A=I2-IBLAT(1)+1 - DSUM=ZEROD - DO 574 I3=IAT(I2),IAT(I2+1)-1 - DSUM=DSUM-AT(I3)*YTMP1(JAT(I3)) - 574 CONTINUE - G2(I2A,I1)=DSUM - 575 CONTINUE - 580 CONTINUE -C /* Now solve g2q2=y for q2 */ -C /* First, input G2 into SPA */ -C -c do 5801 i1=1,iblsz(1,2) -c do 5802 i2=1,iblsz(1,2) -c write(12,*) i1,i2, G2(i1,i2) -c5802 continue -c5801 continue - RCOND(2)=-1.0D0 - CALL DGECO(G2,SIZET,IBLSZ(1,2),KPVT,RCOND(2),Z) -c WRITE(0,*) ' CONDITION ESTIMATES: G1: ',RCOND(1) -c WRITE(0,*) ' CONDITION ESTIMATES: G2: ',RCOND(2) - CALL DGESL(G2,SIZET,IBLSZ(1,2),KPVT,YTMP0) -C -C CALL LUDCMP(G2,IBLSZ(1,2),SIZET,INDX2,DTEMP) -C -C CALL LUBKSB(G2,IBLSZ(1,2),SIZET,INDX2,YTMP0) -C------------------------------------------------------------------------- -c write(8,*) ' IT: ',ITERTN,'G2.............' -c do 582 i=1,iblsz(1,2) -c do 581 j=1,iblsz(1,2) -c write(8,951) itertn, i,j,g2(i,j) -c581 continue -c582 continue -c------------------------------------------------------------------------- */ - -C - DO 605 I1=1,IBLSZ(1,2) - QQ(I1+IBLSZ(1,1))=YTMP0(I1) - 605 CONTINUE -C -C /* Form phat1-At.q2 */ -C /* First, clear ytmp2 */ -C DO 607 I1=1,3000 -C YTMP3(I1)=ZEROD -C607 CONTINUE -C - DO 615 I1=IBLA(1),IBLA(2)-1 - I3=I1-IBLA(1)+1 - YTMP3(I3)=PHAT1(I3) - DO 610 I2=IA(I1),IA(I1+1)-1 - YTMP3(I3)=YTMP3(I3)-A(I2)*QQ(IBLSZ(1,1)+JA(I2)) - 610 CONTINUE - 615 CONTINUE -C /* SOLVE G1q=ytmp3(.) */ - CALL DPOSL(G1,SIZET,IBLSZ(1,1),YTMP3) -C CALL LUBKSB(G1,IBLSZ(1,1),SIZET,INDX,YTMP3) -C - DO 626 I1=1,IBLSZ(1,1) - QQ(I1)=YTMP3(I1) - 626 CONTINUE -c------------------------------------------------------------------------- -C /* NOW qq is complete */ -C .................. SOLVE Sr=Uq .................... -C /* Let dy=dy-r, form r(0) */ - DO 635 I1=IBLAT(1),IBLAT(2)-1 - DSUM=ZEROD - DO 630 I2=IAT(I1),IAT(I1+1)-1 - DSUM=DSUM+AT(I2)*QQ(JAT(I2)) - 630 CONTINUE - I3=I1-IBLAT(1)+1 - IF(I3.GT.IBLSZ(1,2)) WRITE(6,*) 'ERROR2 IN Sr=Uq' - DY(I3)=DY(I3)-QQ(I3+IBLSZ(1,1))-DSUM - 635 CONTINUE - IRWOFS=IBLSZ(1,2) -C /* Form r(i), or equiv */ - DO 660 I1=1,ISCEN -C /* Reload Ch.Fac. of S */ - CALL BLKIO(2,I1,S) - DO 645 I2=IBLAT(I1+1),IBLAT(I1+2)-1 - I2A=I2-IBLAT(I1+1)+1 - DSUM=ZEROD - DO 640 I3=IAT(I2),IAT(I2+1)-1 - DSUM=DSUM+AT(I3)*QQ(JAT(I3)+IACLOF) - 640 CONTINUE - CALL INBI(I2A,DSUM,S) - 645 CONTINUE - CALL SOLVE5(S) SPK06 - DO 650 I2=IRWOFS+1,IRWOFS+IBLSZ(IBLK+1,2) - I3=I2-IRWOFS - DY(I2)=DY(I2)-S(I3) - 650 CONTINUE - IRWOFS=IRWOFS+IBLSZ(I1+1,2) - 660 CONTINUE -C------------------------------------------------------------------------- -c write(8,*) 'ITERTN=',ITERTN -c write(8,903) (i,dy(i),i=1,irwofs+1) -c903 format(2(' dyp(',i3,')=',E13.7)) -c904 format(8(' qq(',i3,')=',F9.5)) -C------------------------------------------------------------------------- - -C /* DONE with dy (EX. PH1) */ - RETURN - END ##### sed > BQ_Decomp/genmat.f <<'#####' 's/^-//' -C -C ************************************************************** -C * * -C * Subroutine BLDMTR: is called only once to initialize * -C * the matrix ADk2At by filling in the * -C * arrays DIAG, AAT. * -C * * -C * called by: Main Program subroutines called: None * -C ************************************************************** -C - SUBROUTINE BLDMTR(iblk,iclofs,IRWOFS,D,ZEROD,EPSILO) - INCLUDE (sizes.h) -C - DOUBLE PRECISION AAT,A,SLACK,DIAG,D(SIZEVEC),LSTPRD,ZEROD - DOUBLE PRECISION ATY(SIZEVEC),B,C,ZCONST,Y,X, EPSILO,SLKBAK -C - INTEGER NROWS,NCOLS,IA,JA,IAAT,JAAT,PTRPRD, LSTDST, I, J -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 /OUTPRO/ LSTPRD(SIZEVEC),PTRPRD(SIZEMAT),LSTDST(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) -C -C To compute TTt, slacks get set to 1 so the computation can proceed. -C Slacks will get overwritten with the correct slacks when the W's -C are processed. /* Find Slacks */ -C - IF(IBLK.EQ.1) GOTO 11 - IF(IBLK.GE.2.AND.IBLK.LE.ISCEN+1) THEN - DO 112 I=ICLOFS+1,ICLOFS+IBLSZ(IBLK,1) - SLACK(I)=1.0D0 - 112 CONTINUE - GOTO 11 - ENDIF -C /* Only form slack for T's once*/ - IF(IBLK.EQ.ISCEN+2) THEN -C /* First, ATY for A */ - DO 3 I1=1,IBLSZ(1,1) - ATY(I1)=ZEROD - I2=I1-IBLA(1)+1 - DO 2 J1=IA(I2),IA(I2+1)-1 - ATY(I1)=ATY(I1)+A(J1)*Y(JA(J1)) - 2 CONTINUE - 3 CONTINUE -C /* Zero out rest of ATY */ - DO 4 I1=IBLSZ(1,1)+1,NCOLS - ATY(I1)=ZEROD - 4 CONTINUE -C /* Now, for T's */ -C /* Shitty code! works- 3-8*/ - DO 7 J2=IACLOF+1,IACLOF+IBLSZ(2,1) - J1=IBLSZ(1,2) - DO 6 I1=2,ISCEN+1 - I2=IBLA(I1)-1+J2-IACLOF -C /*I2 is the row ptr of At */ - DO 5 I3=IA(I2),IA(I2+1)-1 - ATY(J2)=ATY(J2)+A(I3)*Y(J1+JA(I3)) - 5 CONTINUE - J1=J1+IBLSZ(I1,2) - 6 CONTINUE - 7 CONTINUE - DO 8 I1=1,IBLSZ(1,1) - SLACK(I1)=C(I1)-ATY(I1) - 8 CONTINUE - ENDIF -C /* Now compute slack for W's*/ -C /* One for each row of Slack*/ - DO 10 I1=IBLA(IBLK),IBLA(IBLK+1)-1 -C /* J1: element of aty(j1) */ - J1=I1-IBLA(IBLK)+1+ICLOFS - DO 9 I2=IA(I1),IA(I1+1)-1 - ATY(J1)=ATY(J1)+A(I2)*Y(IRWOFS+JA(I2)) - 9 CONTINUE - SLACK(J1)=C(J1)-ATY(J1) -C - 10 CONTINUE -C /* AAT, DIAG zeroed in STJAAT*/ -C /* Below s/b=idiag(iblk+1)-1*/ - 11 SLKBAK=-1.D-01 - DO 110 I1=ICLOFS+1,ICLOFS+IBLSZ(IBLK,1) - IF(DABS(SLACK(I1)).GE.EPSILO) THEN - SLKBAK=-1.0D0*DABS(SLACK(I1)) - GOTO 111 - ENDIF - 110 CONTINUE - 111 CONTINUE - DO 50 I=1, iblsz(iblk,1) - i1=ibla(iblk)-1+i - i2=iclofs+i - IF(DABS(SLACK(I2)).GE.EPSILO) SLKBAK=-1.0D0*DABS(SLACK(I2)) - IF(SLACK(I2).LT.EPSILO.AND.SLACK(I2).GE.-EPSILO) THEN - SLACK(I2)=SLKBAK - ENDIF - D(I2)=1.0D0/slack(i2) - DO 30 J= IA(I1), ia(i1+1)-1 - j1=idiag(iblk)-1+ja(j) - DIAG(j1)=DIAG(j1)+A(J)*A(J)*D(I2)*D(I2) - 30 CONTINUE - DO 40 J=PTRPRD(I), PTRPRD(I+1)-1 -C j1=iblaat(iblk)-1+lstdst(j) - J1=LSTDST(J) - AAT(j1)=AAT(j1)+LSTPRD(J)*D(I2)*D(I2) - 40 CONTINUE - 50 CONTINUE -C - RETURN - END -C ************************************************************ -C * * -C * Subroutine OUTPRD: One computes every nonzero * -C * OUTER PRODUCT A(i,l) x A(j,l), * -C * once in the beginning of the LP procedure. With * -C * these values, this subroutine fills in the arrays * -C * PTRPRD, LSTPRD and LSTDST, which will help at each * -C * iteration, to buils the matrix Bk= ADk2At. * -C * * -C * called by: MAIN program * -C * * -C * subroutines called: None * -C ************************************************************ -C - SUBROUTINE OUTPRD(iblk,error) - INCLUDE (sizes.h) -C - DOUBLE PRECISION LSTPRD,A,AAT,DIAG, LST2 -C - INTEGER PTRPRD,LSTDST,NCOLS,NROWS,IA,JA,IAAT,JAAT - INTEGER L,J,IAUX1,IAUX2,I,K,M,error,iblk,DST2,PTR2 -C - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /CONSTR/ NAMERW(SIZEVEC), RTYPE(SIZEVEC), NAMECL(SIZEVEC) - COMMON /OUTPRO/ LSTPRD(SIZEVEC),PTRPRD(SIZEMAT),LSTDST(SIZEMAT) - COMMON /OUTPR2/ IBLPTR(MAXBLK),IBLLST(MAXBLK),IBLDST(MAXBLK), - 1 LST2(SIZEPRD),DST2(SIZEPRD), PTR2(SIZEVEC) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - 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 PTRPRD: contains pointers to the first nonzero element of -C the outer-product list involving elements in each -C column of A -C LSTPRD: contains the computed outer-products sorted in -C increasing order of the corresponding column index l -C LSTDST: for each entry in LSTPRD, LSTDST contains a pointer -C to the location of the element Bk(i,j) in array AAT -C ---------------------------------------------------------------- - L=1 - DO 130 J=1,iblsz(iblk,1) - j1=ibla(iblk)+J-1 - PTRPRD(J)= L - IAUX2= IA(J1+1) - 2 - IAUX1= IA(J1+1) - 1 -c write(7,*) ' iaux1,2,j1,',iaux1,iaux2,j1,j - IF (IAUX2 .LT. IA(J1)) GOTO 130 - DO 120 I= IA(J1),IAUX2 - DO 110 K=(I+1),IAUX1 - LSTPRD(L)= A(I)*A(K) - j2=iblaat(iblk)+ja(i)-1 - DO 10 M=IAAT(j2),(IAAT(j2+1)-1) -c write(7,*) ' JAAT(',m,')=',JAAT(M),':JA(',K, -c 1 ')=',ja(k) - IF (JAAT(M) .EQ. JA(K)) GOTO 30 - 10 CONTINUE - WRITE(6,20) - 20 FORMAT('ERROR: NO COLUMN IN Bk ASSOCIATES TO OUTER PRD') - WRITE(6,*) ' Trying to find JA(',K,') = ',JA(K),' blk:',iblk - write(6,*) ' col ',namecl(m), ' row : ',namerw(k) - error=1 - GOTO 140 - 30 LSTDST(L)= M - L= L+1 - 110 CONTINUE - 120 CONTINUE - 130 CONTINUE - PTRPRD(iblsz(iblk,1)+1)=L -C - 140 DO 150 I=1,IBLSZ(IBLK,1)+1 - PTR2(IBLPTR(IBLK)+I-1)=PTRPRD(I) - 150 CONTINUE - IBLPTR(IBLK+1)=IBLPTR(IBLK)+IBLSZ(IBLK,1)+1 -C - DO 160 I=1,L-1 - LST2(IBLLST(IBLK)+I-1)=LSTPRD(I) - DST2(IBLDST(IBLK)+I-1)=LSTDST(I) - 160 CONTINUE - IBLLST(IBLK+1)=IBLLST(IBLK)+L-1 - IBLDST(IBLK+1)=IBLDST(IBLK)+L-1 -C - RETURN - END -C -C -C ************************************************************ -C * Subroutine STJAAT: initializes the arrays JAAT AND IAAT * -C * and DIAG,IDIAG * -C * * -C * called by: Main Program subroutines called:COMPCT * -C ************************************************************ -C - SUBROUTINE STJAAT(iblk,ZEROI) - INCLUDE (sizes.h) -C - DOUBLE PRECISION A, AT, aat,diag -C - INTEGER NZTOTL, BEGROW, ENDROW, IAT, JAT, - 1 TOTAUX, NROWS, NCOLS, BEGCOL, ENDCOL, COL, IA, JA, - 1 BEGCL2, ZEROI, JAAT, IAAT, IAUX(4000), - 1 I1, I2, I3, I4 -C - COMMON /ROWISE/ AT(SIZEMAT), IAT(SIZEVEC), JAT(SIZEMAT) - COMMON /INPT1/ A(SIZEMAT), NROWS, NCOLS, IA(SIZEVEC), JA(SIZEMAT) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - COMMON /TMPAAT/ NZTOTL,NIAAT -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 -C -C /* I1=1,NROWS-1 */ - DO 60 I1= IBLAT(IBLK),IBLAT(IBLK+1)-2 - i1a=i1-iblat(iblk)+1 - BEGROW= IAT(I1) - ENDROW= IAT(I1+1) - 1 - TOTAUX= 0 - DO 40 I2= BEGROW, ENDROW - COL= JAT(I2) - BEGCOL= IA(ibla(iblk)+COL-1) - ENDCOL= IA(ibla(iblk)+COL)-1 - IF (I1A .EQ. JA(ENDCOL)) GOTO 40 - DO 20 I3= BEGCOL, ENDCOL - IF (JA(I3) .GT. I1A) THEN - BEGCL2= I3 - DO 10 I4= BEGCL2, ENDCOL -C - TOTAUX= TOTAUX+1 - IAUX(TOTAUX)= JA(I4) - 10 CONTINUE - GOTO 40 - ENDIF - 20 CONTINUE -C - WRITE(6,30) - 30 FORMAT('THIS IS NOT OK. Verify STJAAT subrout.') - 40 CONTINUE -C - CALL COMPCT(IAUX,TOTAUX) -C - DO 50 I2= 1, TOTAUX - NZTOTL= NZTOTL+1 - JAAT(NZTOTL)= IAUX(I2) - AAT(NZTOTL)=0.0D00 - 50 CONTINUE - NIAAT=NIAAT+1 - IAAT(NIAAT)= NZTOTL+1 -C - 60 CONTINUE -C - iblaat(iblk+1)=NIAAT -C /* Initialize IDiag */ - i0=idiag(iblk) - do 70 i=1,iblsz(iblk,2) - diag(i0+i-1)=0.0D00 - 70 continue - idiag(iblk+1)=i0+iblsz(iblk,2) -C - RETURN - END -C -C ************************************************************* -C * * -C * Subroutine SETPH1, introduces a new (constraint in the * -C * primal problem) correponding to adding an * -C * artificial variable in the dual problem. * -C * * -C * called by: Main program * -C * subroutines called: MULTVS, INSERW * -C ************************************************************* -C - SUBROUTINE SETPH1(B,Y,SLKMIN,NROWS,BIGM,BIGMIU) - INCLUDE (sizes.h) -C - DOUBLE PRECISION B(SIZEVEC),Y(SIZEMAT),SLKMIN,ARTFCL,BIGM,BIGMIU -C DOUBLE PRECISION SLACK(SIZEVEC) -C - INTEGER NROWS,I -C -C /* multvs: sc. prod. */ - BIGM=0.0D0 - DO 10 I=1,NROWS - BIGM=BIGM+B(I)*Y(I) - 10 CONTINUE -C - ARTFCL= -2.0D0*SLKMIN - BIGM= DABS(BIGM)/ARTFCL - BIGM= BIGMIU*BIGM -C - Y(NROWS+1)= ARTFCL - B(NROWS+1)= -BIGM -C - RETURN - END -C -C -C *************************************************************** -C * * -C * Subroutine UPDMTR : This subroutine, at each iteration * -C * updates the matrix Bk= ADkDkAt * -C * However the only elements of Dk changing are those * -C * whose relative change in this iteration is greater * -C * then EPSUPD. * -C * * -C * called by: MAIN PROGRAM Subroutines called: None * -C *************************************************************** -C - SUBROUTINE UPDMTR(D,SLACK,EPSUPD) - INCLUDE (sizes.h) -C - DOUBLE PRECISION A,AAT,DIAG,LST2,SLACK(SIZEVEC),D(SIZEMAT) - DOUBLE PRECISION AUX1, AUX2, AUX3, EPSUPD -C - INTEGER NROWS,NCOLS,IA,JA,IAAT,JAAT,PTR2,DST2,JBEG,JEND,I,J - INTEGER JCOL, LTEMP,ICLOFS - character*3 rtype - character*8 NAMERW,NAMECL -C - COMMON /INPT1/ A(SIZEMAT),NROWS,NCOLS,IA(SIZEVEC),JA(SIZEMAT) - COMMON /OUTPR2/ IBLPTR(MAXBLK),IBLLST(MAXBLK),IBLDST(MAXBLK), - 1 LST2(SIZEPRD),DST2(SIZEPRD), PTR2(SIZEVEC) - COMMON /ADK2AT/ AAT(SIZEAAT),DIAG(SIZEVEC),IAAT(SIZEMAT), - 1 JAAT(SIZEAAT) - 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 -C - IBLK=1 - ICLOFS=0 -C - DO 40 I=1,ISCEN+1 -C - IF (I.GT.1) IBLK=ISCEN+I -C - DO 30 J=IBLA(IBLK),IBLA(IBLK+1)-1 - JCOL=J-IBLA(IBLK)+1 - JCOL2=JCOL+ICLOFS -C if(abs(slack(jcol2)).le.1.D-06) then -C write(6,*) ' Error: slack(',jcol2,') = ',slack(jcol2) -C write(6,*) ' Name of ',jcol2,' = ',NAMECL(JCOL2) -C endif - if(abs(D(jcol2)).le.1.D-06) then - write(6,*) ' Error: D(',jcol2,') = ',D(jcol2) - write(6,*) ' Name of ',jcol2,' = ',NAMECL(jcol2) - D(jcol2)=1.0D0 - endif - AUX1=1.0D0/SLACK(JCOL2) - AUX2=DABS(AUX1-D(JCOL2))/DABS(D(JCOL2)) - IF(AUX2.LT.EPSUPD) GOTO 30 - AUX3=(AUX1*AUX1)-(D(JCOL2)*D(JCOL2)) - D(JCOL2)=AUX1 -C - DO 10 K=IA(J),IA(J+1)-1 - JTEMP=IDIAG(IBLK)+JA(K)-1 - DIAG(JTEMP)=DIAG(JTEMP)+A(K)*A(K)*AUX3 - 10 CONTINUE -C - JBEG=PTR2(IBLPTR(IBLK)+JCOL-1) - JEND=PTR2(IBLPTR(IBLK)+JCOL)-1 -C - DO 20 L=JBEG, JEND - LTEMP=DST2(IBLDST(IBLK)-1+L) - AAT(LTEMP)=AAT(LTEMP)+ - 1 LST2(IBLLST(IBLK)-1+L)*AUX3 - 20 CONTINUE -C - 30 CONTINUE - ICLOFS=ICLOFS+IBLSZ(IBLK,1) -C - 40 CONTINUE -C - RETURN - END* -C *************************************************************