/***************************************************************************** * RESOLUTION OF A SET OF NONLINEAR EQUATIONS * * -------------------------------------------------------------------------- * * COPYRIGHT 1991, BY R.S. BAIN * * Originally wrapped by bain@luther.che.wisc.edu on Tue Mar 30 16:27:48 1993 * * * * As of this writing, we do not know how to contact Rod Bain, the * * author of NNES. * * * * -- David M. Gay (dmg@bell-labs.com) * * Bell Labs, Murray Hill * * 8 February 1999 * * * * C++ Release By J-P Moreau, Paris. * * (PART 1/3) * * (www.jpmoreau.fr) * *****************************************************************************/ #include #include //Global variables bool bypass, matsup, wrnsup; int itnum, nfunc; REAL smallb, bigb, smalls, bigs, bigr, nrmpre; //Headers of functions used below REAL Dot_Product(int, REAL *, REAL *); void fcn1(bool *, int, REAL *, REAL *); void fcnevl(bool *, int, int, int, int, REAL, REAL *, REAL *, REAL *); int IMIN(int,int); void innerp (bool, bool *, int, int, int, int, int, int, REAL *, REAL *, REAL *); void Line0(FILE *); void matcop (int, int, int, int, int, int, REAL **, REAL **); REAL MAX(REAL, REAL); REAL MIN(REAL, REAL); void mmulv(int, REAL **, REAL *, REAL *); void Msg(FILE *, char *); void qrsolv(bool, bool *, int, int, int, int, REAL **, REAL *, REAL *, REAL *); void qrupda (bool *, int, int, REAL, REAL **, REAL **, REAL *, REAL *); void rsolv (bool, bool *, int, int, int, int, REAL **, REAL *, REAL *); REAL Sign(REAL, REAL); REAL Sqr(REAL); void twonrm (bool *, int, int, REAL, REAL *, REAL *); void utbmul (int, int, int, int, int, int, REAL **a, REAL **, REAL **); void uvmul (int, int, int, int, REAL **, REAL *, REAL *); void abmul (int nradec, int nraact, int ncbdec, int ncbact, int ncdec, int ncact, REAL **amat, REAL **bmat, REAL **cmat, REAL *arow) { /*----------------------------------------------------------------------- ! FEB. 8, 1991 ! ! MATRIX MULTIPLICATION AB=C ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ! EACH ROW OF MATRIX A IS SAVED AS A COLUMN, AROW, BEfor (E USE. ! ! NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT for ( 1ST INDEX ! NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT for ( 2ND INDEX ! NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ! ! I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ! NCBDEC IS THE NUMBER OF COLUMNS OF B DECLARED ! NCDEC IS THE COMMON DECLARED DIMENSION ! ! MODif (ICATION OF MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ! QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !----------------------------------------------------------------------*/ int i, j, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r; REAL sum; // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=ncact / 32; ncc32r=ncact - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; // REASSIGN ROWS TO VECTOR AROW for (i=1; i<=nraact; i++) { k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; arow[k-31]=amat[i][k-31]; arow[k-30]=amat[i][k-30]; arow[k-29]=amat[i][k-29]; arow[k-28]=amat[i][k-28]; arow[k-27]=amat[i][k-27]; arow[k-26]=amat[i][k-26]; arow[k-25]=amat[i][k-25]; arow[k-24]=amat[i][k-24]; arow[k-23]=amat[i][k-23]; arow[k-22]=amat[i][k-22]; arow[k-21]=amat[i][k-21]; arow[k-20]=amat[i][k-20]; arow[k-19]=amat[i][k-19]; arow[k-18]=amat[i][k-18]; arow[k-17]=amat[i][k-17]; arow[k-16]=amat[i][k-16]; arow[k-15]=amat[i][k-15]; arow[k-14]=amat[i][k-14]; arow[k-13]=amat[i][k-13]; arow[k-12]=amat[i][k-12]; arow[k-11]=amat[i][k-11]; arow[k-10]=amat[i][k-10]; arow[k-9]=amat[i][k-9]; arow[k-8]=amat[i][k-8]; arow[k-7]=amat[i][k-7]; arow[k-6]=amat[i][k-6]; arow[k-5]=amat[i][k-5]; arow[k-4]=amat[i][k-4]; arow[k-3]=amat[i][k-3]; arow[k-2]=amat[i][k-2]; arow[k-1]=amat[i][k-1]; arow[k]=amat[i][k]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; arow[k-15]=amat[i][k-15]; arow[k-14]=amat[i][k-14]; arow[k-13]=amat[i][k-13]; arow[k-12]=amat[i][k-12]; arow[k-11]=amat[i][k-11]; arow[k-10]=amat[i][k-10]; arow[k-9]=amat[i][k-9]; arow[k-8]=amat[i][k-8]; arow[k-7]=amat[i][k-7]; arow[k-6]=amat[i][k-6]; arow[k-5]=amat[i][k-5]; arow[k-4]=amat[i][k-4]; arow[k-3]=amat[i][k-3]; arow[k-2]=amat[i][k-2]; arow[k-1]=amat[i][k-1]; arow[k]=amat[i][k]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; arow[k-7]=amat[i][k-7]; arow[k-6]=amat[i][k-6]; arow[k-5]=amat[i][k-5]; arow[k-4]=amat[i][k-4]; arow[k-3]=amat[i][k-3]; arow[k-2]=amat[i][k-2]; arow[k-1]=amat[i][k-1]; arow[k]=amat[i][k]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; arow[k-3]=amat[i][k-3]; arow[k-2]=amat[i][k-2]; arow[k-1]=amat[i][k-1]; arow[k]=amat[i][k]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; arow[k]=amat[i][k]; } // FIND ENTRY for ( MATRIX C USING COLUMN VECTOR AROW for (j=1; j<=ncbact; j++) { sum=ZERO; k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + arow[k-31]*bmat[k-31][j] + arow[k-30]*bmat[k-30][j] + arow[k-29]*bmat[k-29][j] + arow[k-28]*bmat[k-28][j] + arow[k-27]*bmat[k-27][j] + arow[k-26]*bmat[k-26][j] + arow[k-25]*bmat[k-25][j] + arow[k-24]*bmat[k-24][j]; sum=sum + arow[k-23]*bmat[k-23][j] + arow[k-22]*bmat[k-22][j] + arow[k-21]*bmat[k-21][j] + arow[k-20]*bmat[k-20][j] + arow[k-19]*bmat[k-19][j] + arow[k-18]*bmat[k-18][j] + arow[k-17]*bmat[k-17][j] + arow[k-16]*bmat[k-16][j]; sum=sum + arow[k-15]*bmat[k-15][j] + arow[k-14]*bmat[k-14][j] + arow[k-13]*bmat[k-13][j] + arow[k-12]*bmat[k-12][j] + arow[k-11]*bmat[k-11][j] + arow[k-10]*bmat[k-10][j] + arow[k-9]*bmat[k-9][j] + arow[k-8]*bmat[k-8][j]; sum=sum + arow[k-7]*bmat[k-7][j] + arow[k-6]*bmat[k-6][j] + arow[k-5]*bmat[k-5][j] + arow[k-4]*bmat[k-4][j] + arow[k-3]*bmat[k-3][j] + arow[k-2]*bmat[k-2][j] + arow[k-1]*bmat[k-1][j] + arow[k]*bmat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + arow[k-15]*bmat[k-15][j] + arow[k-14]*bmat[k-14][j] + arow[k-13]*bmat[k-13][j] + arow[k-12]*bmat[k-12][j] + arow[k-11]*bmat[k-11][j] + arow[k-10]*bmat[k-10][j] + arow[k-9]*bmat[k-9][j] + arow[k-8]*bmat[k-8][j]; sum=sum + arow[k-7]*bmat[k-7][j] + arow[k-6]*bmat[k-6][j] + arow[k-5]*bmat[k-5][j] + arow[k-4]*bmat[k-4][j] + arow[k-3]*bmat[k-3][j] + arow[k-2]*bmat[k-2][j] + arow[k-1]*bmat[k-1][j] + arow[k]*bmat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + arow[k-7]*bmat[k-7][j] + arow[k-6]*bmat[k-6][j] + arow[k-5]*bmat[k-5][j] + arow[k-4]*bmat[k-4][j] + arow[k-3]*bmat[k-3][j] + arow[k-2]*bmat[k-2][j] + arow[k-1]*bmat[k-1][j] + arow[k]*bmat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + arow[k-3]*bmat[k-3][j] + arow[k-2]*bmat[k-2][j] + arow[k-1]*bmat[k-1][j] + arow[k]*bmat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + arow[k]*bmat[k][j]; } cmat[i][j]=sum; } //j loop } //i loop } //abmul() void ascalf (int n, REAL epsmch, REAL *fvecc, REAL **jac, REAL *scalef) { /*----------------------------------------------------------------------------- ! FEB. 13, 1991 ! ! THIS SUBROUTINE ESTABLISHES SCALING FACTORS for ( THE RESIDUAL VECTOR ! if ( FUNCTION ADAPTIVE SCALING IS CHOSEN USING INTEGER VARIABLE ITSCLF. ! ! NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE ! UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLICITLY. ! ! SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS OF THE ROWS ! OF THE JACOBIAN AND THE VALUES OF THE CURRENT FUNCTION VECTOR, FVECC. ! ! A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE ! ROOT OF THE MACHINE PRECISION, SQRTEP. !-----------------------------------------------------------------------------*/ REAL amax, sqrtep; int i,j; sqrtep=SQRT(epsmch); // I COUNTS THE ROWS for (i=1; i<=n; i++) { amax=ZERO; // FIND MAXIMUM ENTRY IN ROW I for (j=1; j<=n; j++) amax=MAX(amax, ABS(jac[i][j])); amax=MAX(amax, fvecc[i]); // SET SCALING FACTOR TO A DEFAULT OF ONE if ( ITH ROW IS ZEROS if (amax == ZERO) amax=ONE; amax=MAX(amax, sqrtep); scalef[i]=ONE/amax; } } // ascalf() void ascalx (int n, REAL epsmch, REAL **jac, REAL *scalex) { /*------------------------------------------------------------------------ ! FEB. 13, 1991 ! ! THIS SUBROUTINE ESTABLISHES SCALING FACTORS for ( THE COMPONENET VECTOR ! if ( ADAPTIVE SCALING IS CHOSEN USING INTEGER ITSCLX. ! ! NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE ! UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLICITLY. ! ! SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS ! OF THE COLUMNS OF THE JACOBIAN. ! ! A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE ! ROOT OF THE MACHINE PRECISION, SQRTEP. !-----------------------------------------------------------------------*/ REAL amax, sqrtep; int i,j; sqrtep=SQRT(epsmch); // J COUNTS COLUMNS for (j=1; j<=n; j++) { amax=ZERO; // FIND MAXIMUM ENTRY IN JTH COLUMN for (i=1; i<=n; i++) amax=MAX(amax, ABS(jac[i][j])); // if ( A COLUMN IS ALL ZEROS SET AMAX TO ONE if (amax == ZERO) amax=ONE; scalex[j]=MAX(amax, sqrtep); } } // ascalx() void atamul (int nradec, int ncadec, int nraact, int ncaact, int nrbdec, int ncbdec, REAL **amat, REAL **bmat) { /*---------------------------------------------------------------------- ! FEB. 8, 1991 ! ! MATRIX MULTIPLICATION: A^A=B ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ! ! NRADEC IS NUMBER OF ROWS IN A DECLARED ! NCADEC IS NUMBER OF COLUMNS IN A DECLARED ! NRAACT IS THE LIMIT for ( THE 1ST INDEX IN A ! NCAACT IS THE LIMIT for ( THE 2ND INDEX IN A ! NRBDEC IS NUMBER OF ROWS IN B DECLARED ! NCBDEC IS NUMBER OF COLUMNS IN B DECLARED ! ! MODif (IED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES ! MACKINNON, QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !--------------------------------------------------------------------*/ REAL sum; int i,j,k,kk,ncc4,ncc4r,ncc8,ncc8r,ncc16, ncc16r, ncc32, ncc32r; // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=nraact / 32; ncc32r=nraact - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; // FIND ENTRY IN MATRIX B for (i=1; i<=ncaact; i++) { for (j=i; j<=ncaact; j++) { sum=ZERO; k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[k-31][i]*amat[k-31][j] + amat[k-30][i]*amat[k-30][j] + amat[k-29][i]*amat[k-29][j] + amat[k-28][i]*amat[k-28][j] + amat[k-27][i]*amat[k-27][j] + amat[k-26][i]*amat[k-26][j] + amat[k-25][i]*amat[k-25][j] + amat[k-24][i]*amat[k-24][j]; sum=sum + amat[k-23][i]*amat[k-23][j] + amat[k-22][i]*amat[k-22][j] + amat[k-21][i]*amat[k-21][j] + amat[k-20][i]*amat[k-20][j] + amat[k-19][i]*amat[k-19][j] + amat[k-18][i]*amat[k-18][j] + amat[k-17][i]*amat[k-17][j] + amat[k-16][i]*amat[k-16][j]; sum=sum + amat[k-15][i]*amat[k-15][j] + amat[k-14][i]*amat[k-14][j] + amat[k-13][i]*amat[k-13][j] + amat[k-12][i]*amat[k-12][j] + amat[k-11][i]*amat[k-11][j] + amat[k-10][i]*amat[k-10][j] + amat[k-9][i]*amat[k-9][j] + amat[k-8][i]*amat[k-8][j]; sum=sum + amat[k-7][i]*amat[k-7][j] + amat[k-6][i]*amat[k-6][j] + amat[k-5][i]*amat[k-5][j] + amat[k-4][i]*amat[k-4][j] + amat[k-3][i]*amat[k-3][j] + amat[k-2][i]*amat[k-2][j] + amat[k-1][i]*amat[k-1][j] + amat[k][i]*amat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + amat[k-15][i]*amat[k-15][j] + amat[k-14][i]*amat[k-14][j] + amat[k-13][i]*amat[k-13][j] + amat[k-12][i]*amat[k-12][j] + amat[k-11][i]*amat[k-11][j] + amat[k-10][i]*amat[k-10][j] + amat[k-9][i]*amat[k-9][j] + amat[k-8][i]*amat[k-8][j]; sum=sum + amat[k-7][i]*amat[k-7][j] + amat[k-6][i]*amat[k-6][j] + amat[k-5][i]*amat[k-5][j] + amat[k-4][i]*amat[k-4][j] + amat[k-3][i]*amat[k-3][j] + amat[k-2][i]*amat[k-2][j] + amat[k-1][i]*amat[k-1][j] + amat[k][i]*amat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[k-7][i]*amat[k-7][j] + amat[k-6][i]*amat[k-6][j] + amat[k-5][i]*amat[k-5][j] + amat[k-4][i]*amat[k-4][j] + amat[k-3][i]*amat[k-3][j] + amat[k-2][i]*amat[k-2][j] + amat[k-1][i]*amat[k-1][j] + amat[k][i]*amat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[k-3][i]*amat[k-3][j] + amat[k-2][i]*amat[k-2][j] + amat[k-1][i]*amat[k-1][j] + amat[k][i]*amat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + amat[k][i]*amat[k][j]; } bmat[i][j]=sum; if (i != j) bmat[j][i]=bmat[i][j]; } } } // atamul() void ataov (bool *overfl, int maxexp, int n, int nunit, int output, REAL **a, REAL **b, REAL *scalef) { /*--------------------------------------------------------------------- ! SEPT. 8, 1991 ! ! THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE MATRIX A ! AND MATRIX A. EACH ENTRY IS CHECKED BEfor (E BEING ACCEPTED. ! if ( IT WOULD CAUSE AN OVERFLOW 10^MAXEXP IS INSERTED IN ITS PLACE. !--------------------------------------------------------------------*/ // Labels: e40, e70 REAL eps, sum; int i, j, k; extern FILE *fp; eps=ONE/pow(TEN, maxexp); *overfl=FALSE; for (i=1; i<=n; i++) { for (j=i+1; j<=n; j++) { sum=ZERO; for (k=1; k<=n; k++) { if (log10(ABS(a[k][i])+eps) + log10(ABS(a[k][j])+eps) + TWO*log10(scalef[k]) > maxexp) { *overfl=TRUE; b[i][j]=Sign(pow(TEN,maxexp),a[k][i])*Sign(ONE,a[k][j]); if (output > 2 && (!wrnsup)) { Line0(fp); fprintf(fp," * WARNING: COMPONENT MATRIX-MATRIX PRODUCT SET TO %12.3f *\n", b[i][j]); } goto e40; } sum += a[k][i]*a[k][j]*scalef[k]*scalef[k]; } b[i][j]=sum; b[j][i]=sum; e40:;} // j loop sum=ZERO; for (k=1; k<=n; k++) { if (TWO*(log10(ABS(a[k][i])+eps) + log10(scalef[k])) > maxexp) { *overfl=TRUE; b[i][i]=pow(TEN,maxexp); if (output > 2 && (!wrnsup)) { Line0(fp); fprintf(fp," * WARNING: COMPONENT MATRIX-MATRIX PRODUCT SET TO %12.3f *\n",b[i][i]); } goto e70; } sum += a[k][i]*a[k][i]*scalef[k]*scalef[k]; } b[i][i]=sum; e70:;} // i loop } // ataov() void atbmul (int ncadec, int ncaact, int ncbdec, int ncbact, int ncdec, int ncact, REAL **amat, REAL **bmat, REAL **cmat) { /*--------------------------------------------------------------------- ! FEB. 8, 1991 ! ! MATRIX MULTIPLICATION: A^B=C ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ! ! NCADEC IS 2ND DIM. OF AMAT; NCAACT IS ACTUAL LIMIT for ( 2ND INDEX ! NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT for ( 2ND INDEX ! NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ! ! I.E. NCADEC IS NUMBER OF COLUMNS OF A DECLARED ! NCBDEC IS NUMBER OF COLUMNS OF B DECLARED ! NCDEC IS THE NUMBER OF ROWS IN BOTH A AND B DECLARED ! ! MODIFIED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES ! MACKINNON, QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !---------------------------------------------------------------------*/ int i, j, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r; REAL sum; // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=ncact / 32; ncc32r=ncact - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; // FIND ENTRY IN MATRIX C for (i=1; i<=ncaact; i++) { for (j=1; j<=ncbact; j++) { sum=ZERO; k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[k-31][i]*bmat[k-31][j] + amat[k-30][i]*bmat[k-30][j] + amat[k-29][i]*bmat[k-29][j] + amat[k-28][i]*bmat[k-28][j] + amat[k-27][i]*bmat[k-27][j] + amat[k-26][i]*bmat[k-26][j] + amat[k-25][i]*bmat[k-25][j] + amat[k-24][i]*bmat[k-24][j]; sum=sum + amat[k-23][i]*bmat[k-23][j] + amat[k-22][i]*bmat[k-22][j] + amat[k-21][i]*bmat[k-21][j] + amat[k-20][i]*bmat[k-20][j] + amat[k-19][i]*bmat[k-19][j] + amat[k-18][i]*bmat[k-18][j] + amat[k-17][i]*bmat[k-17][j] + amat[k-16][i]*bmat[k-16][j]; sum=sum + amat[k-15][i]*bmat[k-15][j] + amat[k-14][i]*bmat[k-14][j] + amat[k-13][i]*bmat[k-13][j] + amat[k-12][i]*bmat[k-12][j] + amat[k-11][i]*bmat[k-11][j] + amat[k-10][i]*bmat[k-10][j] + amat[k-9][i]*bmat[k-9][j] + amat[k-8][i]*bmat[k-8][j]; sum=sum + amat[k-7][i]*bmat[k-7][j] + amat[k-6][i]*bmat[k-6][j] + amat[k-5][i]*bmat[k-5][j] + amat[k-4][i]*bmat[k-4][j] + amat[k-3][i]*bmat[k-3][j] + amat[k-2][i]*bmat[k-2][j] + amat[k-1][i]*bmat[k-1][j] + amat[k][i]*bmat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + amat[k-15][i]*bmat[k-15][j] + amat[k-14][i]*bmat[k-14][j] + amat[k-13][i]*bmat[k-13][j] + amat[k-12][i]*bmat[k-12][j] + amat[k-11][i]*bmat[k-11][j] + amat[k-10][i]*bmat[k-10][j] + amat[k-9][i]*bmat[k-9][j] + amat[k-8][i]*bmat[k-8][j]; sum=sum + amat[k-7][i]*bmat[k-7][j] + amat[k-6][i]*bmat[k-6][j] + amat[k-5][i]*bmat[k-5][j] + amat[k-4][i]*bmat[k-4][j] + amat[k-3][i]*bmat[k-3][j] + amat[k-2][i]*bmat[k-2][j] + amat[k-1][i]*bmat[k-1][j] + amat[k][i]*bmat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[k-7][i]*bmat[k-7][j] + amat[k-6][i]*bmat[k-6][j] + amat[k-5][i]*bmat[k-5][j] + amat[k-4][i]*bmat[k-4][j] + amat[k-3][i]*bmat[k-3][j] + amat[k-2][i]*bmat[k-2][j] + amat[k-1][i]*bmat[k-1][j] + amat[k][i]*bmat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[k-3][i]*bmat[k-3][j] + amat[k-2][i]*bmat[k-2][j] + amat[k-1][i]*bmat[k-1][j] + amat[k][i]*bmat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + amat[k][i]*bmat[k][j]; } cmat[i][j]=sum; } } } // atbmul() void atvov (bool *overfl, int maxexp, int n, int nunit, int output, REAL **amat, REAL *bvec, REAL *cvec) { /*--------------------------------------------------------------------- ! FEB. 8 ,1991 ! ! THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE MATRIX A ! AND THE VECTOR B WHERE EACH ENTRY IS CHECKED TO PREVENT OVERFLOWS. !--------------------------------------------------------------------*/ REAL eps, sum; int i, j; char s[80]; extern FILE *fp; eps=ONE/pow(TEN, maxexp); *overfl=FALSE; for (i=1; i<=n; i++) { sum=ZERO; for (j=1; j<=n; j++) { if (log10(ABS(amat[j][i])+eps) + log10(ABS(bvec[j])+eps) > maxexp) { *overfl=TRUE; cvec[i]=Sign(pow(TEN,maxexp), amat[j][i])*Sign(ONE,bvec[j]); if (output > 2 && (!wrnsup)) { Line0(fp); sprintf(s," WARNING: COMPONENT MATRIX-VECTOR PRODUCT SET TO %12.3f", cvec[i]); Msg(fp,s); } return; } } sum += amat[j][i]*bvec[j]; } cvec[i]=sum; } // atvov() void avmul (int nradec, int nraact, int ncdec, int ncact, REAL **amat, REAL *bvec, REAL *cvec) { /*--------------------------------------------------------------------------------- ! FEB. 8, 1991 ! ! MATRIX-VECTOR MULTIPLICATION AB=C ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ! EACH ROW OF MATRIX A IS SAVED AS A COLUMN BEFORE USE. ! ! NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX ! NCDEC IS COMMON DIMENSION OF AMAT & BVEC; NCACT IS ACTUAL LIMIT ! ! I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ! NCDEC IS THE COMMON DECLARED DIMENSION (COLUMNS OF A AND ROWS OF B) ! ! MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ! QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !--------------------------------------------------------------------------------*/ int i, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r; REAL sum; // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=ncact / 32; ncc32r=ncact - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; for (i=1; i<=nraact; i++) { // FIND ENTRY FOR VECTOR C sum=ZERO; k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[i][k-31]*bvec[k-31] + amat[i][k-30]*bvec[k-30] + amat[i][k-29]*bvec[k-29] + amat[i][k-28]*bvec[k-28] + amat[i][k-27]*bvec[k-27] + amat[i][k-26]*bvec[k-26] + amat[i][k-25]*bvec[k-25] + amat[i][k-24]*bvec[k-24]; sum=sum + amat[i][k-23]*bvec[k-23] + amat[i][k-22]*bvec[k-22] + amat[i][k-21]*bvec[k-21] + amat[i][k-20]*bvec[k-20] + amat[i][k-19]*bvec[k-19] + amat[i][k-18]*bvec[k-18] + amat[i][k-17]*bvec[k-17] + amat[i][k-16]*bvec[k-16]; sum=sum + amat[i][k-15]*bvec[k-15] + amat[i][k-14]*bvec[k-14] + amat[i][k-13]*bvec[k-13] + amat[i][k-12]*bvec[k-12] + amat[i][k-11]*bvec[k-11] + amat[i][k-10]*bvec[k-10] + amat[i][k-9]*bvec[k-9] + amat[i][k-8]*bvec[k-8]; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + amat[i][k-15]*bvec[k-15] + amat[i][k-14]*bvec[k-14] + amat[i][k-13]*bvec[k-13] + amat[i][k-12]*bvec[k-12] + amat[i][k-11]*bvec[k-11] + amat[i][k-10]*bvec[k-10] + amat[i][k-9]*bvec[k-9] + amat[i][k-8]*bvec[k-8]; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + amat[i][k]*bvec[k]; } cvec[i]=sum; } } // avmul() void bakdif (bool *overfl, int j, int n, REAL *deltaj, REAL tempj, REAL *fvec, REAL *fvecj1, REAL **jacfdm, REAL *xc) { // FEB. 6, 1991 int i; *deltaj=tempj - xc[j]; fcn1(overfl, n, fvecj1, xc); if (!(*overfl)) for (i=1; i<=n; i++) jacfdm[i][j]=(fvec[i] - fvecj1[i])/(*deltaj); } // bakdif() void bnddif (bool *overfl, int j, int n, REAL epsmch, REAL *boundl, REAL *boundu, REAL *fvecc, REAL *fvecj1, REAL **jacfdm, REAL *xc) { /*------------------------------------------------------------------------------ ! FEB. 15, 1991 ! ! FINITE Dif (FERENCE CALCULATION WHEN THE BOUNDS for ( COMPONENT J ARE SO CLOSE ! THAT NEITHER A for (WARD NOR BACKWARD Dif (FERENCE CAN BE PERfor (MED. !-----------------------------------------------------------------------------*/ int i; REAL eps3q; REAL wv3[10]; eps3q=pow(epsmch,0.75); // STORE CURRENT for (i=1; i<=n; i++) wv3[i]=fvecc[i]; xc[j]=boundu[j]; fcn1 (overfl, n, fvecj1, xc); if (!(*overfl)) { xc[j]=boundl[j]; fcn1 (overfl, n, fvecc, xc); if (!(*overfl)) for (i=1; i<=n; i++) // ENSURE THAT THE JACOBIAN CALCULATION ISN"T JUST NOISE if (fvecj1[i]-fvecc[i] > eps3q) jacfdm[i][j]=(fvecj1[i]-fvecc[i])/(boundu[j]-boundl[j]); else jacfdm[i][j]=ZERO; } for (i=1; i<=n; i++) fvecc[i]=wv3[i]; } // bnddif() void broyfa (bool overch, bool *overfl, bool sclfch, bool sclxch, int maxexp, int n, int nunit, int output, REAL epsmch, REAL **a, REAL *delf, REAL *fvec, REAL *fvecc, REAL **jac, REAL *rdiag, REAL *s, REAL *scalef, REAL *scalex, REAL *t, REAL *w, REAL *xc, REAL *xplus) { /*------------------------------------------------------------------------- ! FEB. 23, 1992 ! ! THE BROYDEN QUASI-NEWTON METHOD IS APPLIED TO THE FACTORED ! FORM OF THE JACOBIAN. ! ! NOTE: T AND W ARE TEMPORARY WORKING VECTORS ONLY. ! ! THE UPDATE OCCURS ONLY IF A SIGNIFICANT CHANGE IN THE JACOBIAN ! WOULD RESULT, I.E., NOT ALL THE VALUES IN VECTOR W ARE LESS THAN ! THE THRESHOLD IN MAGNITUDE. IF THE VECTOR W IS ESSENTIALLY ZERO ) ! THE LOGICAL VARIABLE SKIPUP REMAINS SET AT TRUE. !------------------------------------------------------------------------*/ REAL denom, eps, sum, sqrtep; bool skipup; int i,k; REAL *tmp, **tmp1, **tmp2; void *vmblock = NULL; extern FILE *fp; vmblock = vminit(); tmp1 = (REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); tmp2 = (REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); tmp = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } *overfl=FALSE; eps=ONE/pow(TEN, maxexp); sqrtep=SQRT(epsmch); for (i=1; i<=n; i++) { a[i][i]=rdiag[i]; s[i]=xplus[i] - xc[i]; } // R IS NOW IN THE UPPER TRIANGLE OF A skipup=TRUE; /* THE BROYDEN UPDATE IS CONDENSED INTO THE FOR M A(NEW) = A(OLD) + T S^ THE PRODUCT A*S IS FORMED IN TWO STAGES AS R IS IN THE UPPER TRIANGLE OF MATRIX A AND Q^ IS IN JAC. FIRST MULTIPLY R*S (A IS CONSIDERED UPPER TRIANGULAR) */ uvmul (n, n, n, n, a, s, t); // NOTE: THIS T IS TEMPORARY - IT IS THE T FROM BELOW WHICH // IS SENT TO SUBROUTINE QRUPDA. for (i=1; i<=n; i++) { for (k=1; k<=n; k++) tmp[k]=jac[k][i]; innerp (overch, overfl, maxexp, n, n, n, nunit, output, &sum, tmp, t); w[i]=scalef[i]*(fvec[i]-fvecc[i]) - sum; // TEST TO ENSURE VECTOR W IS NONZERO. ANY VALUE GREATER // THAN THE THRESHOLD WILL SET SKIPUP TO FALSE. if (ABS(w[i]) > sqrtep*scalef[i]*(ABS(fvec[i]) + ABS(fvecc[i]))) skipup=FALSE; else w[i]=ZERO; } // IF W(I)=0 FOR ALL I, THE UPDATE IS SKIPPED. if (!skipup) { // T=Q^W; Q^ IS IN JAC avmul (n, n, n, n, jac, w, t); if (sclxch) for (k=1; k<=n; k++) w[k]=s[k]*scalex[k]; else { for (k=1; k<=n; k++) { tmp1[k][1]=s[k]; tmp2[k][1]=w[k]; } matcop (n, n, 1, 1, n, 1, tmp1, tmp2); for (k=1; k<=n; k++) { s[k]=tmp1[k][1]; w[k]=tmp2[k][1]; } } twonrm (overfl, maxexp, n, epsmch, &denom, w); // IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN if ((*overfl) || log10(denom+eps) > 0.5*maxexp) { if (output > 3) { Line0(fp); Msg(fp," WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF"); Msg(fp," BROYDEN UPDATE."); } goto e10; } else denom=denom*denom; // IF DENOM IS ZERO AVOID DIVIDE BY ZERO AND CONTINUE WITH SAME JACOBIAN if (denom == ZERO) goto e10; // THE SCALED VERSION OF S REPLACES THE ORIGINAL BEFORE BEING SENT TO QRUPDA. for (k=1; k<=n; k++) s[k] *= scalex[k]*scalex[k]/denom; // UPDATE THE QR DECOMPOSITION USING A SERIES OF GIVENS ROTATIONS. qrupda (overfl, maxexp, n, epsmch, a, jac, t, s); // RESET RDIAG AS DIAGONAL OF CURRENT R WHICH IS IN THE UPPER TRIANGE OF A. for (i=1; i<=n; i++) rdiag[i]=a[i][i]; } // UPDATE THE GRADIENT VECTOR, DELF. THE NEW Q^ IS IN JAC. // DELF = (QR)^F = R^Q^F = R^JAC F if (sclfch) for (k=1; k<=n; k++) w[k]=fvec[k]*scalef[k]; else { for (k=1; k<=n; k++) { tmp1[k][1]=fvec[k]; tmp2[k][1]=w[k]; } matcop (n, n, 1, 1, n, 1, tmp1, tmp2); for (k=1; k<=n; k++) { fvec[k]=tmp1[k][1]; w[k]=tmp2[k][1]; } } avmul (n, n, n, n, jac, w, t); for (k=1; k<=n; k++) { tmp1[k][1]=t[k]; tmp2[k][1]=delf[k]; } utbmul (n, n, 1, 1, n, n, a, tmp1, tmp2); for (k=1; k<=n; k++) { t[k]=tmp1[k][1]; delf[k]=tmp2[k][1]; } e10: vmfree(vmblock); } // broyfa() void broyun (bool *overfl, int maxexp, int n, int nunit, int output, REAL epsmch, REAL *fvec, REAL *fvecc, REAL **jac, REAL *scalex, REAL *xc, REAL *xplus) { /*-------------------------------------------------------------------- ! FEB. 23, 1992 ! ! UPDATE THE JACOBIAN USING BROYDEN'S METHOD. !-------------------------------------------------------------------*/ int i, j, k; REAL denom, eps, sqrtep, sum, tempi; REAL *wv1, *tmp1, *tmp2; void *vmblock = NULL; extern FILE *fp; vmblock = vminit(); wv1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); tmp1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); tmp2 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } eps=ONE/pow(TEN, maxexp); sqrtep=SQRT(epsmch); for (k=1; k<=n; k++) wv1[k]=(xplus[k] - xc[k])*scalex[k]; twonrm (overfl, maxexp, n, epsmch, &denom, wv1); // IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN if ((*overfl) || log10(denom+eps) > 0.5*maxexp) { if (output > 3) { Line0(fp); Msg(fp," WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF"); Msg(fp," BROYDEN UPDATE."); } goto e10; } else denom=denom*denom; // IF DENOM IS ZERO, AVOID OVERFLOW, CONTINUE WITH SAME JACOBIAN if (denom == ZERO) goto e10; // UPDATE JACOBIAN BY ROWS for (i=1; i<=n; i++) { for (k=1; k<=n; k++) { tmp1[k]=jac[i][k]; tmp2[k]=xplus[k] - xc[k]; } sum = Dot_Product(n, tmp1, tmp2); tempi=fvec[i] - fvecc[i] - sum; // CHECK TO ENSURE THAT SOME MEANINGFUL CHANGE IS BEING MADE // TO THE APPROXIMATE JACOBIAN; IF NOT, SKIP UPDATING ROW I. if (ABS(tempi) >= sqrtep*(ABS(fvec[i]) + ABS(fvecc[i]))) { tempi=tempi/denom; for (j=1; j<=n; j++) jac[i][j] += tempi*(xplus[j]-xc[j])*scalex[j]*scalex[j]; } } e10: vmfree(vmblock); } // broyun() void cholde (int n, REAL *maxadd, REAL *maxffl, REAL sqrtep, REAL **h, REAL **l) { /*------------------------------------------------------------------------ ! FEB. 23, 1992 ! ! THIS SUBROUTINE FINDS THE CHOLESKY DECOMPOSITION OF THE MATRIX, H, ! AND RETURNS IT IN THE LOWER TRIANGLE OF MATRIX, L. !-----------------------------------------------------------------------*/ REAL minl, minl2, minljj, sum; int i,j,k; REAL *tmp1, *tmp2; void *vmblock = NULL; vmblock = vminit(); tmp1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); tmp2 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } minl=SQRT(sqrtep)*(*maxffl); // MAXFFL EQUALS 0 WHEN THE MATRIX IS KNOWN TO BE POSITIVE DEFINITE if (*maxffl == ZERO) { // FIND SQUARE ROOT OF LARGEST MAGNITUDE DIAGONAL ELEMENT // AND SET MINL2. for (i=1; i<=n; i++) *maxffl=MAX(*maxffl, ABS(h[i][i])); *maxffl=SQRT(*maxffl); minl2=sqrtep*(*maxffl); } // MAXADD CONTAINS THE MAXIMUM AMOUNT WHICH IS IMPLICITLY ADDED // TO ANY DIAGONAL ELEMENT OF MATRIX H. *maxadd=ZERO; for (j=1; j<=n; j++) { for (k=1; k minljj*minljj) // NORMAL CHOLESKY DECOMPOSITION l[j][j]=SQRT(l[j][j]); else { // IMPLICITLY PERTURB DIAGONAL OF H if (minljj < minl2) minljj=minl2; *maxadd=MAX(*maxadd, minljj*minljj-l[j][j]); l[j][j]=minljj; } for (k=j+1; k<=n; k++) l[j+k][j] /= l[j][j]; } vmfree(vmblock); } //cholde() void lsolv (bool, bool *, int, int, int, int, REAL **, REAL *, REAL *); void ltsolv (bool, bool *, int, int, int, int, REAL **, REAL *, REAL *); void chsolv (bool overch, bool *overfl, int maxexp, int n, int nunit, int output, REAL **l, REAL *rhs, REAL *s) { /*-------------------------------------------------------------------- ! ! FEB. 14, 1991 ! ! THIS SUBROUTINE USES for (WARD/BACKWARD SUBSTITUTION TO SOLVE THE ! SYSTEM OF LINEAR EQUATIONS: ! ! (LL^)S=RHS !--------------------------------------------------------------------*/ REAL *wv2; void *vmblock = NULL; vmblock = vminit(); wv2 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } lsolv (overch, overfl, maxexp, n, nunit, output, l, wv2, rhs); ltsolv (overch, overfl, maxexp, n, nunit, output, l, s, wv2); vmfree(vmblock); } // chsolv() void condno (bool overch, bool *overfl, int maxexp, int n, int nunit, int output, REAL *connum, REAL **a, REAL *rdiag) { /*---------------------------------------------------------------------------- ! N.B. Arguments P, PM & Q have been removed. ! ! FEB. 14, 1991 ! ! THIS SUBROUTINE ESTIMATES THE CONDITION NUMBER OF A QR-DECOMPOSED ! MATRIX USING THE METHOD OF CLINE, MOLER, STEWART AND WILKINSON ! (SIAM J. N.A. 16 P368 (1979) ). ! ! if ( A POTENTIAL OVERFLOW IS DETECTED AT ANY POINT ) A CONDITION NUMBER ! EQUIVALENT TO THAT OF A SINGULAR MATRIX IS ASSIGNED BY THE CALLING ! SUBROUTINE. !---------------------------------------------------------------------------*/ int i,j,k; REAL eps, qm, qnorm, qp, temp, tempm; REAL *p, *pm, *q; void *vmblock1 = NULL; vmblock1 = vminit(); p = (REAL *) vmalloc(vmblock1, VEKTOR, n+1, 0); pm = (REAL *) vmalloc(vmblock1, VEKTOR, n+1, 0); q = (REAL *) vmalloc(vmblock1, VEKTOR, n+1, 0); if (! vmcomplete(vmblock1)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } *overfl=FALSE; eps=ONE/pow(TEN, maxexp); *connum=ABS(rdiag[1]); for (j=2; j<=n; j++) { temp=ZERO; for (i=1; i pow(TEN,maxexp-1)) { *overfl=TRUE; goto e10; } temp += ABS(a[i][j]); } temp += ABS(rdiag[j]); *connum=MAX(*connum,temp); } q[1]=ONE/rdiag[1]; for (i=2; i<=n; i++) { if (overch) if (log10(ABS(q[1])+eps) + log10(ABS(a[1][i]) + eps) > maxexp) { *overfl=TRUE; goto e10; } p[i]=a[1][i]*q[1]; } for (j=2; j<=n; j++) { if (overch) if (log10(ABS(p[j])+eps) - log10(ABS(rdiag[j]) + eps) > maxexp) { *overfl=TRUE; goto e10; } qp=(ONE-p[j])/rdiag[j]; qm=(-ONE-p[j])/rdiag[j]; temp=ABS(qp); tempm=ABS(qm); for (i=j+1; i<=n; i++) { if (overch) if (log10(ABS(a[j][i])+eps) + log10(ABS(qm)+eps) > maxexp) { *overfl=TRUE; goto e10; } pm[i]=p[i] + a[j][i]*qm; if (overch) if (log10(ABS(pm[i])+eps) - log10(ABS(rdiag[i])+eps) > maxexp) { *overfl=TRUE; goto e10; } tempm += ABS(pm[i])/ABS(rdiag[i]); if (overch) if (tempm > pow(TEN,maxexp-1)) { *overfl=TRUE; goto e10; } if (overch) if (log10(ABS(a[j][i])+eps) + log10(ABS(qp)+eps) > maxexp) { *overfl=TRUE; goto e10; } p[i] += a[j][i]*qp; if (overch) if (log10(ABS(p[i])+eps) - log10(ABS(rdiag[i])+eps) > maxexp) { *overfl=TRUE; goto e10; } temp += ABS(p[i])/ABS(rdiag[i]); if (overch) if (temp > pow(TEN, maxexp-1)) { *overfl=TRUE; goto e10; } } if (temp >= tempm) q[j]=qp; else { q[j]=qm; for (k=1; k<=n; k++) p[j+k]=pm[j+k]; } } qnorm=ZERO; for (j=1; j<=n; j++) { qnorm += ABS(q[j]); if (overch) if (qnorm > pow(TEN, maxexp-1)) { *overfl=TRUE; goto e10; } } if (log10(*connum) - log10(qnorm) > maxexp) { *overfl=TRUE; goto e10; } *connum = *connum/qnorm; rsolv (overch, overfl, maxexp, n, nunit, output, a, rdiag, q); if (*overfl) goto e10; qnorm=ZERO; for (j=1; j<=n; j++) { qnorm=qnorm + ABS(q[j]); if (overch) if (qnorm > pow(TEN, maxexp-1)) { *overfl=TRUE; goto e10; } } *connum = (*connum)*qnorm; e10: return; } // condno() void delcau (bool cauchy, bool overch, bool *overfl, int itnum, int maxexp, int n, int nunit, int output, REAL *beta, REAL *caulen, REAL *delta, REAL epsmch, REAL *maxstp, REAL *newlen, REAL *sqrtz, REAL **a, REAL *delf, REAL *scalex) { /*------------------------------------------------------------------- ! ! FEB. 23, 1992 ! ! THIS SUBROUTINE ESTABLISHES AN INITIAL TRUST REGION, DELTA, ! if ( ONE IS NOT SPECif (IED BY THE USER AND FINDS THE LENGTH OF ! THE SCALED CAUCHY STEP, CAULEN, AT EACH STEP if ( THE DOUBLE ! DOGLEG OPTION IS BEING USED. ! ! THE USER HAS TWO CHOICES for ( THE INITIAL TRUST REGION: ! 1) BASED ON THE SCALED CAUCHY STEP ! 2) BASED ON THE SCALED NEWTON STEP !------------------------------------------------------------------*/ // Label: e120 #define three 3.0 extern FILE *fp; REAL eps, temp; REAL *wv1; int i,j,k; void *vmblock = NULL; char s[80]; vmblock = vminit(); wv1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } *overfl=FALSE; eps=ONE/pow(TEN, maxexp); /* IF DELTA IS NOT GIVEN EVALUATE IT USING EITHER THE CAUCHY STEP OR THE NEWTON STEP AS SPECIFIED BY THE USER. THE SCALED CAUCHY LENGTH, CAULEN, IS REQUIRED IN TWO CASES. 1) WHEN SELECTED AS THE CRITERION FOR THE INITIAL DELTA 2) IN THE DOUBLE DOGLEG STEP REGARDLESS OF (1) */ if (output > 3) { Line0(fp); Msg(fp," DETERMINATION OF SCALED CAUCHY STEP LENGTH, CAULEN"); } /* FIND FACTOR WHICH GIVES CAUCHY POINT WHEN MULTPLYING STEEPEST DESCENT DIRECTION, DELF. CAULEN = ZETA^1.5/BETA = SQRTZ^3/BETA */ for (k=1; k<=n; k++) wv1[k]=delf[k]/scalex[k]; twonrm (overfl, maxexp, n, epsmch, sqrtz, wv1); if (*overfl) { *caulen=pow(TEN, maxexp); if (output > 4) { Line0(fp); sprintf(s," ZETA SET TO %11.3f TO AVOID OVERFLOW", sqrtz); Msg(fp,s); sprintf(s," SCALED CAUCHY LENGTH, CAULEN SET TO %9.2f TO AVOID OVERFLOW", *caulen); Msg(fp,s); if (itnum = 1) { Line0(fp); Msg(fp," THE PROBLEM SHOULD BE RESCALED OR A NEW STARTING POINT CHOSEN"); Msg(fp," EXECUTION CONTINUES WITH SUBSTITUTIONS AS LISTED ABOVE."); } } } else if (output > 4) { Line0(fp); sprintf(s," SQUARE ROOT OF ZETA, SQRTZ: %12.3f", sqrtz); Msg(fp,s); } // NOTE: THE LOWER TRIANGLE OF MATRIX A NOW CONTAINS THE // TRANSPOSE OF R WHERE A=QR. *beta=ZERO; for (i=1; i<=n; i++) { temp=ZERO; for (j=i; j<=n; j++) { if (overch) if (log10(ABS(a[j][i])+eps) + log10(ABS(delf[j])+eps) > maxexp) { *caulen=SQRT(epsmch); goto e120; } temp += a[j][i]*delf[j]/(scalex[j]*scalex[j]); } *beta += temp*temp; } if (output > 4) { Line0(fp); sprintf(s," BETA: %11.3f NOTE: CAULEN = ZETA^1.5/BETA", *beta); Msg(fp,s); Line0(fp); } // AVOID OVERFLOWS IN FINDING CAULEN if (three*log10(*sqrtz+eps) - log10(*beta+eps) < maxexp && (!(*overfl)) && *beta != ZERO) { // NORMAL DETERMINATION *caulen=*sqrtz*(*sqrtz)*(*sqrtz)/(*beta); /* THIS STEP AVOIDS DIVIDE BY ZERO IN DOGLEG IN THE (RARE) CASE WHERE DELF(I)=0 FOR ALL I BUT THE POINT IS NOT YET A SOLUTION - MOST LIKELY A BAD STARTING ESTIMATE. */ *caulen=MAX(*caulen,ONE/pow(TEN, maxexp)); } else // SUBSTITUTION TO AVOID OVERFLOW *caulen=pow(TEN, maxexp); e120:if (output > 3) { Line0(fp); sprintf(s," SCALED CAUCHY LENGTH, CAULEN: %12.3f", *caulen); Msg(fp,s); } // ESTABLISH INITIAL TRUST REGION IF NEEDED if (*delta < ZERO) { // USE DISTANCE TO CAUCHY POINT OR LENGTH OF NEWTON STEP if (cauchy) *delta = MIN(*caulen, *maxstp); else *delta = MIN(*newlen, *maxstp); if (output > 3) { Line0(fp); sprintf(s," INITIAL TRUST REGION SIZE, DELTA: %12.3f", *delta); Msg(fp,s); } } vmfree(vmblock); } // delcau() void deufls (bool *abort, bool deuflh, bool *geoms, bool overch, bool *overfl, bool *qnfail, bool *qrsing, bool *restrt, bool sclfch, bool *sclxch, int *acpcod, int *acptcr, int *contyp, int *itnum, int *jupdm, int maxexp, int maxlin, int n, int *nfunc, int nunit, int output, int qnupdm, int *stopcr, REAL alpha, REAL confac, REAL delfts, REAL epsmch, REAL fcnmax, REAL *fcnnew, REAL fcnold, REAL *lambda, REAL newmax, REAL *sbrnrm, REAL sigma, REAL **a, REAL **astore, REAL *boundl, REAL *boundu, REAL *delf, REAL *fvec, REAL *hhpi, REAL **jac, REAL *rdiag, REAL *rhs, REAL *s, REAL *sbar, REAL *scalef, REAL *scalex, REAL *sn, REAL *xc, REAL *xplus) { /*--------------------------------------------------------------------------- ! FEB. 23, 1992 ! ! THIS SUBROUTINE CONDUCTS A LINE SEARCH IN THE NEWTON DIRECTION ! IF NO CONSTRAINTS ARE VIOLATED. IF THE FIRST TRIAL IS A FAILURE ! THERE ARE TWO TYPES OF LINE SEARCH AVAILABLE: ! 1) REDUCE THE RELAXATION FACTOR, LAMBDA, TO SIGMA*LAMBDA ! WHERE SIGMA IS USER-SPECif (IED (GEOMETRIC Line0 SEARCH) ! 2) AT THE FIRST STEP MINIMIZE A QUADRATIC THROUGH THE OBJECTIVE ! FUNCTION AT THE CURRENT POINT AND TRIAL ESTIMATE (WHICH MUST BE A ! FAILURE) WITH INITIAL SLOPE DELFTS. AT SUBSEQUENT STEPS MINIMIZE ! A CUBIC THROUGH THE OBJECTIVE FUNCTION AT THE TWO MOST RECENT ! FAILURES AND THE CURRENT POINT, AGAIN USING THE INITIAL SLOPE, ! DELFTS. ! ! CONVIO INDICATES A CONSTRAINT VIOLATION BY ONE OR MORE COMPONENTS ! FRSTST INDICATES THE FIRST STEP IN THE LINE SEARCH. ! ! RATIO RATIO OF PROPOSED STEP LENGTH IN (I)TH DIRECTION ! TO DISTANCE FROM (I)TH COMPONENT TO BOUNDARY VIOLATED ! RATIOM MINIMUM OF RATIOS FOR ALL CONSTAINTS VIOLATED !--------------------------------------------------------------------------*/ #define point1 0.1 #define point5 0.5 #define three 3.0 REAL acubic, bcubic, disc, dlftsm, eps, factor, fplpre, lampre, lamtmp, lambda2, ratio, ratiom; REAL *wv2; int i, ii, j, k; bool convio, frstst; void *vmblock = NULL; char ss[80]; extern FILE *fp; vmblock = vminit(); wv2 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } frstst=TRUE; *overfl=FALSE; *abort = FALSE; eps=ONE/pow(TEN, maxexp); *lambda=ONE; //added by JPM for (k=1; k<=maxlin; k++) { ratiom=ONE; convio=FALSE; // FIND TRIAL POINT AND CHECK IF CONSTRAINTS VIOLATED // (IF CONTYP IS NOT EQUAL TO ZERO). for (i=1; i<=n; i++) { /* NOTE: WV2 MARKS VIOLATIONS. WV2(I) CHANGES TO 1 FOR LOWER BOUND VIOLATIONS AND TO 2 FOR UPPER BOUND VIOLATIONS. CONSTRAINT VIOLATIONS CAN ONLY OCCUR AT THE FIRST STEP. */ wv2[i]=-ONE; xplus[i]=xc[i] + (*lambda)*sn[i]; if (contyp > 0 && (frstst)) if (xplus[i] < boundl[i]) { convio=TRUE; wv2[i]=ONE; } else if (xplus[i] > boundu[i]) { convio=TRUE; wv2[i]=TWO; } else wv2[i]=-ONE; } /* IF CONSTRAINTS ARE VIOLATED FIRST REDUCE THE STEP SIZES FOR THE VIOLATING COMPONENTS TO OBTAIN A FEASIBLE POINT. IF THE DIRECTION TO THIS MODIFIED POINT IS NOT A DESCENT DIRECTION OR IF THE MODIFIED STEP DOES NOT LEAD TO AN ACCEPTABLE POINT THEN RETURN TO THE NEWTON DIRECTION AND START A LINE SEARCH AT A FEASIBLE POINT WHERE THE COMPONENT WHICH HAS THE SMALLEST VALUE OF RATIO (DEFINED BELOW) IS TAKEN TO "CONFAC" OF THE DISTANCE TO THE BOUNDARY. DEFAULT VALUE OF CONFAC IS 0.95. */ if (convio) { if (output > 3) { Line0(fp); sprintf(ss," LINE SEARCH STEP: %3d", k); Msg(fp,ss); Line0(fp); sprintf(ss," LAMBDA FOR ATTEMPTED STEP: %12.3f", *lambda); Msg(fp,ss); Msg(fp," CONSTRAINT VIOLATED"); Line0(fp); Msg(fp," TRIAL ESTIMATES (VIOLATIONS MARKED)."); Line0(fp); for (i=1; i<=n; i++) if (wv2[i] > ZERO) { sprintf(ss," XPLUS(%3d) = %12.3f *", i, xplus[i]); Msg(fp,ss); } else { sprintf(ss," XPLUS(%3d) = %12.3f", i, xplus[i]); Msg(fp,ss); } } for (i=1; i<=n; i++) if (wv2[i] > ZERO) { // FIND RATIO for ( THIS VIOLATING COMPONENT } if (wv2[i] == ONE) ratio=-(xc[i]-boundl[i])/(xplus[i]-xc[i]); else if (wv2[i] == TWO) ratio=(boundu[i]-xc[i])/(xplus[i]-xc[i]); // NOTE: THIS Line0 IS for ( OUTPUT PURPOSES ONLY } wv2[i]=ratio; ratiom=MIN(ratiom,ratio); if (ratio > point5) s[i]=confac*ratio*(*lambda)*sn[i]; else // WITHIN BUFFER ZONE - ONLY TAKE 1/2 // OF THE STEP YOU WOULD TAKE OTHERWISE. s[i]=point5*confac*ratio*(*lambda)*sn[i]; // ESTABLISH MODIFIED TRIAL POINT. xplus[i]=xc[i] + s[i]; } else { // FOR NONVIOLATORS XPLUS REMAINS UNCHANGED BUT THE COMPONENT // OF S IS LOADED TO CHECK THE DIRECTIONAL DERIVATIVE. s[i]=(*lambda)*sn[i]; } if (output > 3) { Line0(fp); Msg(fp," NEW S AND XPLUS VECTORS (WITH RATIOS for ( VIOLATIONS)"); Line0(fp); Msg(fp," NOTE: RATIOS ARE RATIO OF LENGTH TO BOUNDARY FROM CURRENT"); Msg(fp," X VECTOR TO MAGNITUDE OF CORRESPONDING PROPOSED STEP."); Line0(fp); for (i=1; i<=n; i++) if (wv2[i] < ZERO) { sprintf(ss," S(%3d) = %12.3f XPLUS(%3d) = %12.3f", i,s[i],i,xplus[i]); Msg(fp,ss); } else { sprintf(ss," S(%3d) = %12.3f XPLUS(%3d) = %12.3f %11.3f",i,s[i],i,xplus[i],wv2[i]); Msg(fp,ss); } Line0(fp); sprintf(ss," MINIMUM OF RATIOS, RATIOM: %12.3f", ratiom); Msg(fp,ss); } // CHECK DIRECTIONAL DERIVATIVE FOR MODIFIED POINT, DLFTSM innerp (overch, overfl, maxexp, n, n, n, nunit, output, &dlftsm, delf, s); *overfl=FALSE; if (output > 3) { Line0(fp); sprintf(ss," INNER PRODUCT OF DELF AND S FOR MODIFIED S: %12.3f", dlftsm); Msg(fp,ss); } // IF INNER PRODUCT IS POSITIVE RETURN TO NEWTON DIRECTION // AND CONDUCT A LINE SEARCH WITHIN THE FEASIBLE REGION. if (dlftsm > ZERO) { if (output > 3) { Msg(fp," DELFTS > 0 START LINE SEARCH AT LAMBDA = CONFAC*LAMBDA*RATIOM"); Msg(fp," NOTE: NO TRIAL POINT WAS EVALUATED AT THIS STEP OF LINE SEARCH."); } // THE STARTING POINT IS SET AT JUST INSIDE THE MOST VIOLATED BOUNDARY. *lambda=confac*ratiom*(*lambda); // LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE } } // if (convio) // NO CONSTRAINTS VIOLATED - EVALUATE RESIDUAL VECTOR AT NEW POINT fcn1 (overfl, n, fvec, xplus); *nfunc = *nfunc + 1; // CHECK FOR OVERFLOW IN FUNCTION VECTOR EVALUATION. // IF SO, REDUCE STEP LENGTH AND CONTINUE LINE SEARCH. if (*overfl) { if (output > 3) Msg(fp," OVERFLOW IN FUNCTION VECTOR - STEP LENGTH REDUCED."); // FORCE STEP TO BE WITHIN CONSTRAINTS - DON'T CALL // THIS THE FIRST STEP, I.E. FRSTST STAYS AT TRUE. if (convio) *lambda=ratiom*confac*(*lambda); else *lambda=sigma*(*lambda); // LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE } // EVALUATE (POSSIBLY SCALED) OBJECTIVE FUNCTION AT NEW POINT fcnevl (overfl, maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef); if (*overfl) { if (output > 3) Msg(fp," OVERFLOW IN OBJECTIVE FUNCTION - STEP LENGTH REDUCED."); // FORCE STEP TO BE WITHIN CONSTRAINTS - DON"T CALL // THIS THE FIRST STEP, I.E. FRSTST STAYS AT TRUE. if (convio) *lambda=ratiom*confac*(*lambda); else *lambda=sigma*(*lambda); } /* IF DEUFLHARD"S METHOD IS BEING USED FOR EITHER RELAXATION FACTOR INITIALIZATION OR THE SECOND ACCEPTANCE CRITERION, EVALUATE SBAR. EVALUATION METHOD DEPENDS UPON WHETHER THE JACOBIAN WAS PERTURBED IN THE SOLUTION OF THE LINEAR SYSTEM. LOGICAL VARIABLE QRSING IS TRUE IF PERTURBATION TOOK PLACE. */ if ((deuflh) || *acptcr == 12) { if (*qrsing) { /* FOR M -J^F AS RIGHT HAND SIDE - METHOD DEPENDS ON WHETHER QNUPDM EQUALS 0 OR 1 IF A QUASI-NEWTON UPDATE IS BEING USED. IF JUPDM IS 0 THEN THE NEWTON STEP HAS BEEN FOUND IN SUBROUTINE NSTPUN. */ if (jupdm == 0 || qnupdm == 0) { // UNSCALED JACOBIAN IN MATRIX JAC for (i=1; i<=n; i++) if (sclfch) wv2[i]=-fvec[i]*scalef[i]*scalef[i]; else wv2[i]=-fvec[i]; mmulv (n, jac, wv2, rhs); } else { // R IN UPPER TRIANGLE OF A PLUS RDIAG AND Q^ IN JAC // - FROM QR DECOMPOSITION OF SCALED JACOBIAN. for (i=1; i<=n; i++) { wv2[i]=ZERO; for (j=1; j<=n; j++) wv2[i] -= jac[i][j]*fvec[j]*scalef[j]; } rhs[1]=rdiag[1]*wv2[1]; for (j=2; j<=n; j++) { rhs[j]=ZERO; for (i=1; i 3) { Line0(fp); sprintf(ss," LINE SEARCH STEP: %3d", k); Msg(fp,ss); Line0(fp); sprintf(ss," LAMBDA FOR ATTEMPTED STEP: %12.3f",*lambda); Msg(fp,ss); Line0(fp); if (!convio) Msg(fp," NEW COMPONENT/FCN VECTORS (XPLUS(I) = XC(I) + LAMBDA*SN(I))"); else Msg(fp," NEW FUNCTION VECTORS AT MODIFIED POINT"); Line0(fp); for (i=1; i<=n; i++) { if (ABS(fvec[i])<100000) sprintf(ss," XPLUS(%3d) = %12.3f FVEC(%3d) = %12.3f",i, xplus[i], i, fvec[i]); else sprintf(ss," XPLUS(%3d) = %12.3f FVEC(%3d) = %e",i, xplus[i], i, fvec[i]); Msg(fp,ss); } Line0(fp); if (!sclfch) { sprintf(ss," OBJECTIVE FUNCTION VALUE AT XPLUS: %12.3f", fcnnew); Msg(fp,ss); } else { sprintf(ss," SCALED OBJECTIVE FUNCTION VALUE AT XPLUS: %12.3f", fcnnew); Msg(fp,ss); } if (ABS(fcnmax+alpha*(*lambda)*delfts) < 100000) sprintf(ss," FCNMAX + ALPHA*LAMBDA*DELFTS: %12.3f", fcnmax+alpha*(*lambda)*delfts); else sprintf(ss," FCNMAX + ALPHA*LAMBDA*DELFTS: %e", fcnmax+alpha*(*lambda)*delfts); Msg(fp,ss); if ((deuflh) || *acptcr == 12) if (itnum > 0) if (!(*sclxch)) { Line0(fp); Msg(fp," DEUFLHARD SBAR VECTOR:"); Line0(fp); for (i=1; i<=n; i++) { sprintf(ss," SBAR(%3d) = %e", i, sbar[i]); Msg(fp,ss); } } else { Line0(fp); Msg(fp," DEUFLHARD SBAR VECTOR IN SCALE DX UNITS:"); Line0(fp); for (i=1; i<=n; i++) { sprintf(ss," SBAR(%3d) = %e", i, scalex[i]*sbar[i]); Msg(fp,ss); } } if (*acptcr == 12) { Line0(fp); if (!(*sclxch)) { sprintf(ss," VALUE OF SBRNRM AT XPLUS: ..................%12.3f", sbrnrm); Msg(fp,ss); } else { sprintf(ss," VALUE OF SCALED SBRNRM AT XPLUS: ...........%12.3f", sbrnrm); Msg(fp,ss); } sprintf(ss," NEWMAX:: ...................................%12.3f", newmax); Msg(fp,ss); } } // *** CHECK FOR ACCEPTABLE STEP *** if (*fcnnew < fcnmax + alpha*(*lambda)*delfts) { *acpcod=1; // NOTE: STEP WILL BE ACCEPTED REGARDLESS OF NEXT TEST. // THIS SECTION IS FOR BOOKKEEPING ONLY. if (*acptcr == 12) if (*sbrnrm < newmax) *acpcod=12; vmfree(vmblock); return; } if (*acptcr == 12 && *sbrnrm < newmax) { *acpcod=2; vmfree(vmblock); return; } // FAILURE OF STEP ACCEPTANCE TEST if (convio) *lambda=confac*ratiom*(*lambda); // LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE if (*lambda == ZERO) { if (output > 0) { Line0(fp); Msg(fp," LAMBDA IS 0.0: NO PROGRESS POSSIBLE - CHECK BOUNDS OR START."); } *abort=TRUE; vmfree(vmblock); return; } if (*geoms) // GEOMETRIC LINE SEARCH *lambda=sigma*(*lambda); else { if (frstst) { frstst=FALSE; // FIND MINIMUM OF QUADRATIC AT FIRST STEP lambda2 = *lambda*(*lambda); lamtmp=-lambda2*delfts/(TWO*(*fcnnew-fcnold-(*lambda)*delfts)); if (output > 4) { Line0(fp); sprintf(ss," TEMPORARY LAMBDA FROM QUADRATIC MODEL: %11.3f", lamtmp); Msg(fp,ss); } } else { // FIND MINIMUM OF CUBIC AT SUBSEQUENT STEPS factor=ONE/(*lambda-lampre); lambda2 = *lambda*(*lambda); if (lambda2 == ZERO ) *lambda=sigma*(*lambda); // NOTE: IF THIS LAMBDA^2 WAS ZERO ANY SUBSEQUENT // LAMBDA^2 WILL ALSO BE ZERO. acubic=factor*((ONE/lambda2)*(*fcnnew-fcnold-*lambda*delfts)-((ONE/lampre*(*lambda))*(fplpre-fcnold-lampre*delfts))); bcubic=factor*((-lampre/lambda2)*(*fcnnew-fcnold-*lambda*delfts)+((*lambda/lampre*(*lambda)* (fplpre-fcnold-lampre*delfts)))); if (TWO*log10(ABS(bcubic)+eps) > maxexp) { lamtmp=sigma*(*lambda); if (output > 4) { Line0(fp); Msg(fp," POTENTIAL OVERFLOW IN CALCULATING TRIAL LAMBDA FROM"); Msg(fp," CUBIC MODEL - LAMBDA SET TO SIGMA*LAMBDA."); } } else { disc=bcubic*bcubic - three*acubic*delfts; if (ABS(acubic) <= epsmch) lamtmp=-delfts/(TWO*bcubic); else if (disc < ZERO) lamtmp=sigma*(*lambda); else lamtmp=(-bcubic + SQRT(disc))/(three*acubic); if (output > 4) { Line0(fp); sprintf(ss," TEMPORARY LAMBDA FROM CUBIC MODEL : .....%11.3f", lamtmp); } if (lamtmp > sigma*(*lambda)) { lamtmp=sigma*(*lambda); if (output > 4) { Line0(fp); Msg(fp," LAMTMP TOO LARGE - REDUCED TO SIGMA*LAMBDA."); } } } // 3rd else } //2nd else } // 1st else lampre=*lambda; fplpre=*fcnnew; if (lamtmp < point1*(*lambda)) { if (output > 4) { Line0(fp); Msg(fp," LAMTMP TOO SMALL - INCREASED TO 0.1*LAMBDA."); } *lambda=point1*(*lambda); } else { if (output > 4 && lamtmp != sigma*(*lambda)) { Line0(fp); Msg(fp," LAMTMP WITHIN LIMITS - LAMBDA SET TO LAMTMP."); } *lambda=lamtmp; } } // k loop! // FAILURE IN LINE SEARCH *acpcod=0; /* IF A QUASI-NEWTON STEP HAS FAILED IN THE LINE SEARCH, SET QNFAIL TO TRUE ANS RETURN TO SUBROUTINE NNES. THIS WILL CAUSE THE JACOBIAN TO BE RE-EVALUATED EXPLICITLY AND A LINE SEARCH IN THE NEW DIRECTION CONDUCTED. */ if (!restrt) { *qnfail=TRUE; vmfree(vmblock); return; } // FALL THROUGH MAIN LOOP - WARNING GIVEN if (output > 2 && (!wrnsup)) { Line0(fp); sprintf(ss," WARNING: %3d CYCLES COMPLETED IN LINE SEARCH WITHOUT SUCCESS", maxlin); Msg(fp,ss); } if (*stopcr == 2 || *stopcr == 3) { *stopcr=12; Line0(fp); Msg(fp," STOPPING CRITERION RESET FROM 2 TO 12 TO AVOID HANGING."); } vmfree(vmblock); } // deufls() void fcnevl (bool *overfl, int maxexp, int n, int nunit, int output, REAL epsmch, REAL *fcnnew, REAL *fvec, REAL *scalef) { /*--------------------------------------------------------------------- ! FEB. 23, 1992 ! ! THE OBJECTIVE FUNCTION, FCNNEW, DEFINED BY: ! ! FCNNEW = 1/2(SCALEF*FVEC^SCALEF*FVEC) ! ! IS EVALUATED BY THIS SUBROUTINE. !--------------------------------------------------------------------*/ REAL eps; int k; REAL wv1[10]; char s[10]; extern FILE *fp; *overfl=FALSE; eps=ONE/pow(TEN, maxexp); for (k=1; k<=n; k++) wv1[k]=fvec[k]*scalef[k]; twonrm (overfl, maxexp, n, epsmch, fcnnew, wv1); // IF AN OVERFLOW WOULD OCCUR SUBSTITUTE A LARGE VALUE FOR FCNNEW if ((*overfl) || TWO*log10(*fcnnew+eps) > maxexp) { *overfl=TRUE; *fcnnew=pow(TEN, maxexp); if (output > 2 && (!wrnsup)) { Line0(fp); sprintf(s," WARNING: TO AVOID OVERFLOW, OBJECTIVE FUNCTION SET TO: %11.3f", *fcnnew); Msg(fp,s); } return; } *fcnnew = *fcnnew*(*fcnnew)/TWO; } //fcnevl() void fordif (bool *overfl, int j, int n, REAL deltaj, REAL *fvec, REAL *fvecj1, REAL **jacfdm, REAL *xc) { /*--------------------------------------------- ! FEB. 6, 1991 !--------------------------------------------*/ int k; fcn1 (overfl, n, fvecj1, xc); if (!(*overfl)) for (k=1; k<=n; k++) jacfdm[k][j]=(fvecj1[k] - fvec[k])/deltaj; } void innerp (bool overch, bool *overfl, int maxexp, int ldima, int ldimb, int n, int nunit, int output, REAL *dtpro, REAL *a, REAL *b) { /*------------------------------------------------------------------- ! FEB. 14, 1991 ! ! THIS SUBROUTINE FINDS THE INNER PRODUCT OF TWO VECTORS, A AND B. ! IF OVERCH IS FALSE, UNROLLED LOOPS ARE USED. ! ! LDIMA IS THE DIMENSION OF A ! LDIMB IS THE DIMENSION OF B ! N IS THE DEPTH INTO A AND B THE INNER PRODUCT IS DESIRED. ! (USUALLY LDIMA=LDIMB=N) !------------------------------------------------------------------*/ REAL eps; int i, k, kk, ng4, ng4r, ng8, ng8r, ng16, ng16r, ng32, ng32r; char s[80]; extern FILE *fp; eps=ONE/pow(TEN, maxexp); *overfl=FALSE; *dtpro=ZERO; if (overch) for (i=1; i<=n; i++) { if (log10(ABS(a[i])+eps) + log10(ABS(b[i])+eps) > maxexp) { *overfl=TRUE; *dtpro=Sign(pow(TEN,maxexp),a[i])*Sign(ONE,b[i]); if (output > 2 && (!wrnsup)) { Line0(fp); sprintf(s," * WARNING: TO AVOID OVERFLOW, INNER PRODUCT SET TO %12.3f", *dtpro); Msg(fp,s); } return; } *dtpro = *dtpro + a[i]*b[i]; } else { // SET NUMBER OF GROUPS OF EACH SIZE ng32=n / 32; ng32r=n - 32*ng32; ng16=ng32r / 16; ng16r=ng32r - 16*ng16; ng8=ng16r / 8; ng8r=ng16r - 8*ng8; ng4=ng8r / 4; ng4r=ng8r - 4*ng4; // FIND INNER PRODUCT k=0; if (ng32 > 0) for (kk=1; kk<=ng32; kk++) { k=k+32; *dtpro=*dtpro + a[k-31]*b[k-31] + a[k-30]*b[k-30] + a[k-29]*b[k-29] + a[k-28]*b[k-28] + a[k-27]*b[k-27] + a[k-26]*b[k-26] + a[k-25]*b[k-25] + a[k-24]*b[k-24]; *dtpro=*dtpro + a[k-23]*b[k-23] + a[k-22]*b[k-22] + a[k-21]*b[k-21] + a[k-20]*b[k-20] + a[k-19]*b[k-19] + a[k-18]*b[k-18] + a[k-17]*b[k-17] + a[k-16]*b[k-16]; *dtpro=*dtpro + a[k-15]*b[k-15] + a[k-14]*b[k-14] + a[k-13]*b[k-13] + a[k-12]*b[k-12] + a[k-11]*b[k-11] + a[k-10]*b[k-10] + a[k-9]*b[k-9] + a[k-8]*b[k-8]; *dtpro=*dtpro + a[k-7]*b[k-7] + a[k-6]*b[k-6] + a[k-5]*b[k-5] + a[k-4] *b[k-4] + a[k-3]*b[k-3] + a[k-2]*b[k-2] + a[k-1]*b[k-1] + a[k]*b[k]; } if (ng16 > 0) for (kk=1; kk<=ng16; kk++) { k=k+16; *dtpro=*dtpro + a[k-15]*b[k-15] + a[k-14]*b[k-14] + a[k-13]*b[k-13] + a[k-12]*b[k-12] + a[k-11]*b[k-11] + a[k-10]*b[k-10] + a[k-9]*b[k-9] + a[k-8]*b[k-8]; *dtpro=*dtpro + a[k-7]*b[k-7] + a[k-6]*b[k-6] + a[k-5]*b[k-5] + a[k-4]*b[k-4] + a[k-3]*b[k-3] + a[k-2]*b[k-2] + a[k-1]*b[k-1] + a[k]*b[k]; } if (ng8 > 0) for (kk=1; kk<=ng8; kk++) { k=k+8; *dtpro=*dtpro + a[k-7]*b[k-7] + a[k-6]*b[k-6] + a[k-5]*b[k-5] + a[k-4]*b[k-4] + a[k-3]*b[k-3] + a[k-2]*b[k-2] + a[k-1]*b[k-1] + a[k]*b[k]; } if (ng4 > 0) for (kk=1; kk<=ng4; kk++) { k=k+4; *dtpro=*dtpro + a[k-3]*b[k-3] + a[k-2]*b[k-2] + a[k-1]*b[k-1] + a[k]*b[k]; } if (ng4r > 0) for (kk=1; kk<=ng4r; kk++) { k++; *dtpro=*dtpro + a[k]*b[k]; } } } // innerp() void jacrot (bool *overfl, int i, int maxexp, int n, REAL arot, REAL brot, REAL epsmch, REAL **a, REAL **jac) { /*--------------------------------------------------------- ! FEB. 11, 1991 ! ! JACOBI (OR GIVENS) ROTATION. !--------------------------------------------------------*/ REAL c,denom,s,w,y; REAL hold[3]; //index 0 not used int j, ldhold; if (arot == ZERO) { c=ZERO; s=-Sign(ONE, brot); } else { hold[1]=arot; hold[2]=brot; ldhold=2; twonrm (overfl, maxexp, ldhold, epsmch, &denom, hold); c=arot/denom; s=-brot/denom; } for (j=i; j<=n; j++) { y=a[i][j]; w=a[i+1][j]; a[i][j]=c*y - s*w; a[i+1][j]=s*y + c*w; } for (j=1; j<=n; j++) { y=jac[i][j]; w=jac[i+1][j]; jac[i][j]=c*y - s*w; jac[i+1][j]=s*y + c*w; } } // jacrot() void lsolv (bool overch, bool *overfl, int maxexp, int n, int nunit, int output, REAL **l, REAL *b, REAL *rhs) { /*---------------------------------------------------------------------------- ! FEB. 14, 1991 ! ! THIS SUBROUTINE SOLVES: ! ! LB=RHS ! ! WHERE L IS TAKEN FROM THE CHOLESKY DECOMPOSITION ! RHS IS A GIVEN RIGHT HAND SIDE WHICH IS NOT OVERWRITTEN ! B IS THE SOLUTION VECTOR ! ! FRSTER IS USED for ( OUTPUT PURPOSES ONLY. !---------------------------------------------------------------------------*/ //Label: e30 int i, j, jstar; REAL eps, maxlog, sum, tmplog; bool frster; char s[80]; extern FILE *fp; frster=TRUE; *overfl=FALSE; eps=ONE/pow(TEN,maxexp); if (overch) if (log10(ABS(rhs[1])+eps) - log10(ABS(l[1][1])+eps) > maxexp) { *overfl=TRUE; b[1]=Sign(pow(TEN,maxexp), rhs[1])*Sign(ONE,l[1][1]); if (output > 3) { Line0(fp); sprintf(s," WARNING: COMPONENT 1 SET TO %12.3f", b[1]); Msg(fp,s); } goto e30; } b[1]=rhs[1]/l[1][1]; e30: for (i=2; i<=n; i++) { if (overch) { // CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW. maxlog=log10(ABS(rhs[i])+eps) - log10(ABS(l[i][i])+eps); jstar=0; for (j=1; j maxlog) { jstar=j; maxlog=tmplog; } } // IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE TERM WITH CORRECT SIGN. if (maxlog > maxexp) { *overfl=TRUE; if (jstar == 0) b[i]=Sign(pow(TEN,maxexp), rhs[i])*Sign(ONE,l[i][i]); else b[i]=-Sign(pow(TEN,maxexp), l[i][jstar])*Sign(ONE,b[jstar])*Sign(ONE,l[i][i]); if (frster) { frster=FALSE; Line0(fp); } if (output > 3) { sprintf(s," WARNING: COMPONENT %3d SET TO %12.3f", i, b[1]); Msg(fp,s); } } } // SUM FOR EACH TERM, ORDERING OPERATIONS TO MINIMIZE POSSIBILITY OF OVERFLOW. sum=ZERO; for (j=1; j maxexp) { *overfl=TRUE; y[n]=Sign(pow(TEN,maxexp),b[n])*Sign(ONE,l[n][n]); if (output > 3) { frster=FALSE; Line0(fp); sprintf(s," WARNING: COMPONENT %3d SET TO %11.3f", n, y[n]); Msg(fp,s); } goto e30; } y[n]=b[n]/l[n][n]; e30: for (i=n-1; i>0; i--) { if (overch) { // CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW maxlog=log10(ABS(b[i])+eps) - log10(ABS(l[i][i])+eps); jstar=0; for (j=i+1; j<=n; j++) { tmplog=log10(ABS(l[j][i])+eps) + log10(ABS(y[j])+eps) - log10(ABS(l[i][i])+eps); if (tmplog > maxlog) { jstar=j; maxlog=tmplog; } } // IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE TERM WITH CORRECT SIGN. if (maxlog > maxexp) { *overfl=TRUE; if (jstar == 0) y[i]=Sign(pow(TEN,maxexp),b[i])*Sign(ONE,l[i][i]); else y[i]=-Sign(pow(TEN,maxexp),l[jstar][i])*Sign(ONE,y[jstar])*Sign(ONE,l[i][i]); if (frster) { frster=FALSE; Line0(fp); } if (output > 3) { sprintf(s," WARNING: COMPONENT %3d SET TO %11.3f", i, y[i]); Msg(fp,s); } } } // SUM FOR EACH TERM ORDERING OPERATIONS TO MINIMIZE POSSIBILITY OF OVERFLOW. sum=ZERO; for (j=i+1; j<=n; j++) sum += (MIN(ABS(l[j][i]),ABS(y[j]))/l[i][i])*(MAX(ABS(l[j][i]),ABS(y[j])))*Sign(ONE,l[j][i])*Sign(ONE,y[j]); y[i]=b[i]/l[i][i] - sum; } } // ltsolv() void matcop (int nradec, int nraact, int ncadec, int ncaact, int nrbdec, int ncbdec, REAL **amat, REAL **bmat) { /*----------------------------------------------------------------- ! SEPT. 15, 1991 ! ! COPY A CONTINGUOUS RECTANGULAR PORTION OF ONE MATRIX ! INTO ANOTHER (ELEMENT [1][1] MUST BE INCLUDED). ! ! NRADEC IS 1ST DIMENSION OF AMAT, NRAACT IS LIMIT OF 1ST INDEX ! NCADEC IS 2ND DIMENSION OF AMAT, NCAACT IS LIMIT OF 2ND INDEX ! NRBDEC IS 1ST DIMENSION OF BMAT ! NCBDEC IS 2ND DIMENSION OF BMAT !----------------------------------------------------------------*/ int j, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r; // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=nraact / 32; ncc32r=nraact - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; for (j=1; j<=ncaact; j++) { // COPY ENTRIES INTO MATRIX B BY COLUMN k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; bmat[k-31][j]=amat[k-31][j]; bmat[k-30][j]=amat[k-30][j]; bmat[k-29][j]=amat[k-29][j]; bmat[k-28][j]=amat[k-28][j]; bmat[k-27][j]=amat[k-27][j]; bmat[k-26][j]=amat[k-26][j]; bmat[k-25][j]=amat[k-25][j]; bmat[k-24][j]=amat[k-24][j]; bmat[k-23][j]=amat[k-23][j]; bmat[k-22][j]=amat[k-22][j]; bmat[k-21][j]=amat[k-21][j]; bmat[k-20][j]=amat[k-20][j]; bmat[k-19][j]=amat[k-19][j]; bmat[k-18][j]=amat[k-18][j]; bmat[k-17][j]=amat[k-17][j]; bmat[k-16][j]=amat[k-16][j]; bmat[k-15][j]=amat[k-15][j]; bmat[k-14][j]=amat[k-14][j]; bmat[k-13][j]=amat[k-13][j]; bmat[k-12][j]=amat[k-12][j]; bmat[k-11][j]=amat[k-11][j]; bmat[k-10][j]=amat[k-10][j]; bmat[k-9][j]=amat[k-9][j]; bmat[k-8][j]=amat[k-8][j]; bmat[k-7][j]=amat[k-7][j]; bmat[k-6][j]=amat[k-6][j]; bmat[k-5][j]=amat[k-5][j]; bmat[k-4][j]=amat[k-4][j]; bmat[k-3][j]=amat[k-3][j]; bmat[k-2][j]=amat[k-2][j]; bmat[k-1][j]=amat[k-1][j]; bmat[k][j]=amat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; bmat[k-15][j]=amat[k-15][j]; bmat[k-14][j]=amat[k-14][j]; bmat[k-13][j]=amat[k-13][j]; bmat[k-12][j]=amat[k-12][j]; bmat[k-11][j]=amat[k-11][j]; bmat[k-10][j]=amat[k-10][j]; bmat[k-9][j]=amat[k-9][j]; bmat[k-8][j]=amat[k-8][j]; bmat[k-7][j]=amat[k-7][j]; bmat[k-6][j]=amat[k-6][j]; bmat[k-5][j]=amat[k-5][j]; bmat[k-4][j]=amat[k-4][j]; bmat[k-3][j]=amat[k-3][j]; bmat[k-2][j]=amat[k-2][j]; bmat[k-1][j]=amat[k-1][j]; bmat[k][j]=amat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; bmat[k-7][j]=amat[k-7][j]; bmat[k-6][j]=amat[k-6][j]; bmat[k-5][j]=amat[k-5][j]; bmat[k-4][j]=amat[k-4][j]; bmat[k-3][j]=amat[k-3][j]; bmat[k-2][j]=amat[k-2][j]; bmat[k-1][j]=amat[k-1][j]; bmat[k][j]=amat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; bmat[k-3][j]=amat[k-3][j]; bmat[k-2][j]=amat[k-2][j]; bmat[k-1][j]=amat[k-1][j]; bmat[k][j]=amat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; bmat[k][j]=amat[k][j]; } } } // matcop() void qrsolv(bool overch, bool *overfl, int maxexp, int n, int nunit, int output, REAL **a, REAL *hhpi, REAL *rdiag, REAL *b) { /*-------------------------------------------------------------------------- ! FEB. 2, 1991 ! ! THIS SUBROUTINE SOLVES ! ! [QR]X=B ! ! WHERE Q AND R ARE OBTAINED FROM THE QR DECOMPOSITION ! B IS A GIVEN RIGHT HAND SIDE WHICH IS OVERWRITTEN ! ! R IS CONTAINED IN THE STRICT UPPER TRIANGLE OF ! MATRIX A AND THE VECTOR RDIAG ! Q IS "CONTAINED" IN THE LOWER TRIANGLE OF MATRIX A ! ! FRSTOV INDICATES FIRST OVERFLOW - USED ONLY TO SET BORDER for ( OUTPUT !-------------------------------------------------------------------------*/ bool frstov; REAL eps, tau; int i,j; char s[80]; extern FILE *fp; eps=ONE/pow(TEN, maxexp); frstov=TRUE; *overfl=FALSE; // MULTIPLY RIGHT HAND SIDE BY Q THEN SOLVE USING R STORED IN MATRIX A for (j=1; j maxexp) { *overfl=TRUE; tau=Sign(pow(TEN,maxexp), a[i][j])*Sign(ONE, b[j]); } tau += a[i][j]*b[i]/hhpi[j]; } for (i=j; i<=n; i++) { if (overch) if (log10(ABS(tau)+eps) + log10(ABS(a[i][j])+eps) > maxexp) { *overfl=TRUE; b[i]=-Sign(pow(TEN,maxexp),tau)*Sign(ONE,a[i][j]); if (output > 2 && (!wrnsup)) { if (frstov) Line0(fp); sprintf(s," WARNING: COMPONENT %3d SET TO %11.3f IN QRSOLV BEFORE RSOLV", i, b[i]); Msg(fp,s); } } b[i] -= tau*a[i][j]; } } rsolv (overch, overfl, maxexp, n, nunit, output, a, rdiag, b); } // qrsolv() void qrupda (bool *overfl, int maxexp, int n, REAL epsmch, REAL **a, REAL **jac, REAL *u, REAL *v) { /*---------------------------------------------------------------- ! FEB. 12, 1991 ! ! UPDATE QR DECOMPOSITION USING A SERIES OF GIVENS ROTATIONS. !---------------------------------------------------------------*/ REAL eucnrm; REAL hold[3]; //index 0 not used int i, k, l, ldhold; // REPLACE SUBDIAGONAL WITH ZEROS SO THAT WHEN R IS MULTIPLIED BY GIVENS // (JACOBI) ROTATIONS THE SUBDIAGONAL ELEMENTS DO NOT AFFECT THE OUTCOME for (i=2; i<=n; i++) a[i][i-1]=ZERO; // FIND LARGEST K FOR WHICH U[K] DOES NOT EQUAL ZERO k=n; for (l=1; l<=n; l++) if (u[k] == ZERO) if (k > 1) k--; else return; else return; /* MULTIPLY UV^ BY A SERIES OF ROTATIONS SO THAT ALL BUT THE TOP ROW IS MADE ZERO (THEORETICALLY THIS IS WHAT HAPPENS ALTHOUGH THIS MATRIX ISN'T ACTUALLY FORMED). */ for (i=k-1; i>0; i--) { jacrot (overfl, i, maxexp, n, u[i], u[i+1], epsmch, a, jac); if (u[i] == ZERO) // THIS STEP JUST AVOIDS ADDING ZERO u[i]=ABS(u[i+1]); else { hold[1]=u[i]; hold[2]=u[i+1]; ldhold=2; twonrm (overfl, maxexp, ldhold, epsmch, &eucnrm, hold); u[i]=eucnrm; } } // ADD THE TOP ROW TO THE TOP ROW OF A - THIS FORMS THE // UPPER HESSENBERG MATRIX. for (i=1; i<=n; i++) a[1][i] += u[1]*v[i]; // FORM THE UPPER TRIANGULAR R MATRIX BY A SERIES OF ROTATIONS // TO ZERO OUT THE SUBDIAGONALS. for (i=1; i maxexp) { *overfl=TRUE; b[n]=Sign(pow(TEN,maxexp),b[n])*Sign(ONE,rdiag[n]); if (output > 2 && (!wrnsup)) { frstov=FALSE; Line0(fp); sprintf(s," WARNING: COMPONENT %3d SET TO %12.3f", n, b[n]); Msg(fp,s); } goto e30; } b[n] /= rdiag[n]; e30:for (i=n-1; i>0; i--) { if (overch) { // CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW. maxlog=log10(ABS(b[i])+eps) - log10(ABS(rdiag[i])+eps); jstar=0; for (j=i+1; j<=n; j++) { tmplog=log10(ABS(a[i][j])+eps) + log10(ABS(b[j])+eps) - log10(ABS(rdiag[i])+eps); if (tmplog > maxlog) { jstar=j; maxlog=tmplog; } } // IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE TERM WITH CORRECT SIGN. if (maxlog > maxexp) { *overfl=TRUE; if (jstar == 0) b[i]=Sign(pow(TEN,maxexp), b[i])*Sign(ONE, rdiag[i]); else b[i]=-Sign(pow(TEN,maxexp), a[i][jstar])*Sign(ONE, b[jstar])*Sign(ONE, rdiag[i]); if (frstov) { frstov=FALSE; Line0(fp); } if (output > 2 && (!wrnsup)) { sprintf(s," WARNING: COMPONENT %3d SET TO %12.3f", i, b[i]); Msg(fp,s); } } } // SUM FOR EACH TERM ORDERING OPERATIONS TO MINIMIZE POSSIBILITY OF OVERFLOW. sum=ZERO; for (j=i+1; j<=n; j++) sum += (MIN(ABS(a[i][j]),ABS(b[j]))/rdiag[i])*(MAX(ABS(a[i][j]),ABS(b[j])))*Sign(ONE,a[i][j])*Sign(ONE,b[j]); b[i]=b[i]/rdiag[i]-sum; } } // rsolv() void setup (bool *absnew, bool *cauchy, bool *deuflh, bool *geoms, bool *linesr, bool *newton, bool *overch, int * acptcr, int *itsclf, int *itsclx, int *jactyp, int *jupdm, int *maxexp, int *maxit, int *maxns, int *maxqns, int *minqns, int n, int *narmij, int *niejev, int *njacch, int *output, int *qnupdm, int *stopcr, int *supprs, int *trupdm, REAL *alpha, REAL *confac, REAL *delta, REAL *delfac, REAL *epsmch, REAL *etafac, REAL *fdtolj, REAL *ftol, REAL *lam0, REAL *mstpf, REAL *nsttol, REAL *omega, REAL *ratiof, REAL *sigma, REAL *stptol, REAL *boundl, REAL *boundu, REAL *scalef, REAL *scalex, char *help) { /*----------------------------------------------------------------------------- ! DEC. 7, 1991 ! ! FUNCTION SETUP ASSIGNS DEFAULT VALUES TO ALL REQUISITE PARAMETERS. !----------------------------------------------------------------------------*/ REAL temp, xmax; int i, ibeta, it, maxebb, minebb; extern nfetot; // BOOLEAN VALUES *absnew=FALSE; bypass=FALSE; *cauchy=FALSE; *deuflh=TRUE; *geoms= TRUE; *linesr=TRUE; matsup=FALSE; *newton=FALSE; *overch=FALSE; wrnsup=FALSE; // INTEGER VALUES *acptcr=12; *itsclf=0; *itsclx=0; *jactyp=1; *jupdm=0; *maxit=100; //max. number of iterations *maxns=50; *maxqns=10; *minqns=7; *narmij=1; nfetot=0; *niejev=1; *njacch=1; *output=2; // set to 4 to have more outputs *qnupdm=1; *stopcr=12; *supprs=0; *trupdm=0; // REAL VALUES *alpha=1E-04; *confac=0.95; *delta=-1.0; *delfac=2.0; *etafac=0.2; *lam0=1.0; *mstpf=1E3; *omega=0.1; *ratiof=0.7; *sigma=0.5; // CHARACTER VARIABLE strcpy(help,"NONE"); /* NOTE: NOTATIONAL CHANGES IN CALLING PROGRAM FROM MACHAR 1) EPSMCH DENOTES MACHINE EPSILON 2) MINEBB DENOTES MINIMUM EXPONENT BASE BETA 3) MAXEBB DENOTES MAXIMUM EXPONENT BASE BETA */ // Values extracted from Fortran 90 version it = 53; ibeta = 2; minebb = -1021; maxebb = 1024; *epsmch = 2.220446049250313E-16; xmax = 1.7976931348E+300; *maxexp = 308; // VALUES FOR TWO-NORM CALCULATIONS smallb=pow(ONE*ibeta,(minebb+1) / 2); bigb=pow(ONE*ibeta,(maxebb-it+1) / 2); smalls=pow(ONE*ibeta,(minebb-1) / 2); bigs=pow(ONE*ibeta,(maxebb+it-1) / 2); bigr=xmax; // SET STOPPING CRITERIA PARAMETERS *fdtolj=1E-06; *ftol=pow(*epsmch,0.333); *nsttol=(*ftol)*(*ftol); *stptol=*nsttol; // VECTOR VALUES temp=-pow(TEN, *maxexp); for (i=1; i<=n; i++) { boundl[i]=temp; boundu[i]=-temp; scalef[i]=ONE; scalex[i]=ONE; } } // setup() void twonrm (bool *overfl, int maxexp, int n, REAL epsmch, REAL *eucnrm, REAL *v) { /*---------------------------------------------------------------------- ! FEB. 23 ,1992 ! ! THIS SUBROUTINE EVALUATES THE EUCLIDEAN NORM OF A VECTOR, V. ! IT FOLLOWS THE ALGORITHM OF J.L. BLUE IN ACM TOMS V4 15 (1978) ! BUT USES SLIGHTLY Dif (FERENT CUTS. THE CONSTANTS IN COMMON BLOCK ! NNES_5 ARE CALCULATED IN THE SUBROUTINE MACHAR OR ARE PROVIDED ! BY THE USER IN THE DRIVER. !---------------------------------------------------------------------*/ REAL abig, absvi, amed, asmall, sqrtep, ymax, ymin; int i; *overfl=FALSE; sqrtep=SQRT(epsmch); asmall=ZERO; amed=ZERO; abig=ZERO; for (i=1; i<=n; i++) { // ACCUMULATE SUMS OF SQUARES IN ONE OF THREE ACCULULATORS, // ABIG, AMED AND ASMALL, DEPENDING ON THEIR SIZES. absvi=ABS(v[i]); // THIS COMPARISON RESTRICTS THE MAXIMUM VALUE OF AMED TO BE // B/N => CANNOT SUM SO THAT AMED OVERFLOWS. if (absvi > bigb/(ONE*n)) // THIS /ISOR OF 10N RESTRICTS ABIG FROM (PATHALOGICALLY) // OVERFLOWING FROM SUMMATION. abig += v[i]/Sqr(TEN*n*bigs); else if (absvi < smalls) asmall += Sqr(v[i]/smalls); else amed += v[i]*v[i]; } if (abig > ZERO) { // IF OVERFLOW WOULD OCCUR ASSIGN BIGR AS NORM AND SIGNAL TO // CALLING SUBROUTINE VIA OVERFL. if (log10(abig)/TWO + log10(bigs) + ONE + log10(ONE*n) > maxexp) { *eucnrm=bigr; *overfl=TRUE; return; } // IF AMED IS POSITIVE IT COULD CONTRIBUTE TO THE NORM - // DETERMINATION IS DELAYED UNTIL LATER TO SAVE REPEATING CODE. if (amed > ZERO) { ymin=MIN(SQRT(amed),TEN*n*bigs*SQRT(abig)); ymax=MAX(SQRT(amed),TEN*n*bigs*SQRT(abig)); } else { // AMED DOESN'T CONTRIBUTE AND ASMALL WON'T MATTER IF // ABIG IS NONZERO - FIND NORM USING ABIG AND RETURN. *eucnrm=TEN*n*bigs*SQRT(abig); return; } } else if (asmall > ZERO) if (amed > ZERO) { ymin=MIN(SQRT(amed),smalls*SQRT(asmall)); ymax=MAX(SQRT(amed),smalls*SQRT(asmall)); } else { // ABIG AND AMED ARE ZERO SO USE ASMALL ONLY. *eucnrm=smalls*SQRT(asmall); return; } else { *eucnrm=SQRT(amed); return; } if (ymin < sqrtep*ymax) // SMALLER PORTION DOES NOT CONTRIBUTE TO NORM. *eucnrm=ymax; else // SMALLER PORTION CONTRIBUTES TO NORM. *eucnrm=ymax*SQRT((ONE+ymin/ymax)*(ONE+ymin/ymax)); } // twonrm() void update (int mnew, int *mold, int n, int *trmcod, REAL *fcnnew, REAL *fcnold, REAL *fvec, REAL *fvecc, REAL *xc, REAL *xplus) { /*----------------------------------------------------------------------- ! FEB. 9, 1991 ! ! THIS SUBROUTINE RESETS CURRENT ESTIMATES OF SOLUTION AND UPDATES THE ! OBJECTIVE FUNCTION VALUE, M (USED TO SET HOW MANY PREVIOUS VALUES TO ! LOOK AT IN THE NON- MONOTONIC COMPARISONS) AND THE TERMINATION CODE, ! TRMCOD. !----------------------------------------------------------------------*/ int i; *fcnold=*fcnnew; *mold=mnew; *trmcod=0; for (i=1; i<=n; i++) { fvecc[i]=fvec[i]; xc[i]=xplus[i]; } } // update() void utbmul (int ncadec, int ncaact, int ncbdec, int ncbact, int ncdec, int ncact, REAL **amat, REAL **bmat, REAL **cmat) { /*-------------------------------------------------------------------------- ! FEB. 8, 1991 ! ! MATRIX MULTIPLICATION: A^B=C WHERE A IS UPPER TRIANGULAR ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ! ! NCADEC IS 2ND DIM. OF AMAT; NCAACT IS ACTUAL LIMIT for ( 2ND INDEX ! NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT for ( 2ND INDEX ! NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ! ! I.E. NCADEC IS NUMBER OF COLUMNS OF A DECLARED ! NCBDEC IS NUMBER OF COLUMNS OF B DECLARED ! NCDEC IS THE NUMBER OF ROWS IN BOTH A AND B DECLARED ! ! MODif (ICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ! QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !-------------------------------------------------------------------------*/ int i, j, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r, nend; REAL sum; for (i=1; i<=ncaact; i++) { // FIND NUMBER OF GROUPS OF SIZE 32, 16... nend=IMIN(i,ncact); // THIS ADJUSTMENT IS REQUIRED WHEN NCACT IS LESS THAN NCAACT ncc32=nend / 32; ncc32r=nend - 32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; // FIND ENTRY IN MATRIX C for (j=1; j<=ncbact; j++) { sum=ZERO; k=0; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[k-31][i]*bmat[k-31][j]+amat[k-30][i]*bmat[k-30][ j]+amat[k-29][i]*bmat[k-29][j]+amat[k-28][i]*bmat[k-28][j]+ amat[k-27][i]*bmat[k-27][j]+amat[k-26][i]*bmat[k-26][j]+ amat[k-25][i]*bmat[k-25][j]+amat[k-24][i]*bmat[k-24][j]; sum=sum + amat[k-23][i]*bmat[k-23][j]+amat[k-22][i]*bmat[k-22][ j]+amat[k-21][i]*bmat[k-21][j]+amat[k-20][i]*bmat[k-20][j]+ amat[k-19][i]*bmat[k-19][j]+amat[k-18][i]*bmat[k-18][j]+ amat[k-17][i]*bmat[k-17][j]+amat[k-16][i]*bmat[k-16][j]; sum=sum + amat[k-15][i]*bmat[k-15][j]+amat[k-14][i]*bmat[k-14][ j]+amat[k-13][i]*bmat[k-13][j]+amat[k-12][i]*bmat[k-12][j]+ amat[k-11][i]*bmat[k-11][j]+amat[k-10][i]*bmat[k-10][j]+ amat[k-9][i]*bmat[k-9][j]+amat[k-8][i]*bmat[k-8][j]; sum=sum + amat[k-7][i]*bmat[k-7][j]+amat[k-6][i]*bmat[k-6][j]+ amat[k-5][i]*bmat[k-5][j]+amat[k-4][i]*bmat[k-4][j]+amat[k-3][ i]*bmat[k-3][j]+amat[k-2][i]*bmat[k-2][j]+amat[k-1][i]* bmat[k-1][j]+amat[k][i]*bmat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + amat[k-15][i]*bmat[k-15][j]+amat[k-14][i]*bmat[k-14][ j]+amat[k-13][i]*bmat[k-13][j]+amat[k-12][i]*bmat[k-12][j]+ amat[k-11][i]*bmat[k-11][j]+amat[k-10][i]*bmat[k-10][j]+ amat[k-9][i]*bmat[k-9][j]+amat[k-8][i]*bmat[k-8][j]; sum=sum + amat[k-7][i]*bmat[k-7][j]+amat[k-6][i]*bmat[k-6][j]+ amat[k-5][i]*bmat[k-5][j]+amat[k-4][i]*bmat[k-4][j]+amat[k-3][ i]*bmat[k-3][j]+amat[k-2][i]*bmat[k-2][j]+amat[k-1][i]* bmat[k-1][j]+amat[k][i]*bmat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[k-7][i]*bmat[k-7][j]+amat[k-6][i]*bmat[k-6][j]+ amat[k-5][i]*bmat[k-5][j]+amat[k-4][i]*bmat[k-4][j]+amat[k-3][ i]*bmat[k-3][j]+amat[k-2][i]*bmat[k-2][j]+amat[k-1][i]* bmat[k-1][j]+amat[k][i]*bmat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[k-3][i]*bmat[k-3][j]+amat[k-2][i]*bmat[k-2][j]+ amat[k-1][i]*bmat[k-1][j]+amat[k][i]*bmat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + amat[k][i]*bmat[k][j]; } cmat[i][j]=sum; } } } // utbmul() void uvmul (int nradec, int nraact, int ncdec, int ncact, REAL **amat, REAL *bvec, REAL *cvec) { /*------------------------------------------------------------------------------ ! FEB. 8, 1991 ! ! MATRIX-VECTOR MULTIPLICATION: AB=C WHERE A IS UPPER TRIANGULAR ! ! VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ! EACH ROW OF MATRIX A IS SAVED AS A COLUMN BEFORE USE. ! ! NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX ! NCDEC IS COMMON DIMENSION OF AMAT & BVEC; NCACT IS ACTUAL LIMIT ! ! I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ! NCDEC IS THE COMMON DECLARED DIMENSION (COLUMNS OF A AND ROWS OF B) ! ! MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ! QUEEN"S UNIVERSITY, KINGSTON, ONTARIO, CANADA !------------------------------------------------------------------------------*/ int i, k, kk, ncc4, ncc4r, ncc8, ncc8r, ncc16, ncc16r, ncc32, ncc32r; REAL sum; for (i=1; i<=nraact; i++) { // FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32=(ncact - (i-1)) / 32; ncc32r=(ncact - (i-1))-32*ncc32; ncc16=ncc32r / 16; ncc16r=ncc32r - 16*ncc16; ncc8=ncc16r / 8; ncc8r=ncc16r - 8*ncc8; ncc4=ncc8r / 4; ncc4r=ncc8r - 4*ncc4; // FIND ENTRY FOR VECTOR C sum=ZERO; k=i-1; if (ncc32 > 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[i][k-31]*bvec[k-31] + amat[i][k-30]*bvec[k-30] + amat[i][k-29]*bvec[k-29] + amat[i][k-28]*bvec[k-28] + amat[i][k-27]*bvec[k-27] + amat[i][k-26]*bvec[k-26] + amat[i][k-25]*bvec[k-25] + amat[i][k-24]*bvec[k-24]; sum=sum + amat[i][k-23]*bvec[k-23] + amat[i][k-22]*bvec[k-22] + amat[i][k-21]*bvec[k-21] + amat[i][k-20]*bvec[k-20] + amat[i][k-19]*bvec[k-19] + amat[i][k-18]*bvec[k-18] + amat[i][k-17]*bvec[k-17] + amat[i][k-16]*bvec[k-16]; sum=sum + amat[i][k-15]*bvec[k-15] + amat[i][k-14]*bvec[k-14] + amat[i][k-13]*bvec[k-13] + amat[i][k-12]*bvec[k-12] + amat[i][k-11]*bvec[k-11] + amat[i][k-10]*bvec[k-10] + amat[i][k-9]*bvec[k-9] + amat[i][k-8]*bvec[k-8]; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; kk++) { k=k+16; sum=sum + amat[i][k-15]*bvec[k-15] + amat[i][k-14]*bvec[k-14] + amat[i][k-13]*bvec[k-13] + amat[i][k-12]*bvec[k-12] + amat[i][k-11]*bvec[k-11] + amat[i][k-10]*bvec[k-10] + amat[i][k-9]*bvec[k-9] + amat[i][k-8]*bvec[k-8]; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[i][k-7]*bvec[k-7] + amat[i][k-6]*bvec[k-6] + amat[i][k-5]*bvec[k-5] + amat[i][k-4]*bvec[k-4] + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[i][k-3]*bvec[k-3] + amat[i][k-2]*bvec[k-2] + amat[i][k-1]*bvec[k-1] + amat[i][k]*bvec[k]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum=sum + amat[i][k]*bvec[k]; } cvec[i]=sum; } // i loop } // uvmul() //end of file unnes.cpp