C ALGORITHM 741, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 21, NO. 1, MARCH, 1995, P. 20-25. C c*** bbqr.f * Least Squares Solution of a * Linear, Bordered, Block-diagonal System of Equations * * This package of 6 subroutines solves by least squares an overdeter- * mined, full-rank, single-bordered, block-diagonal system of equations. * The system is solved in a sequential fashion, using a series of * orthogonal transformations in the form of Householder reflections. * For this documentation to make sense, the user must understand how * the terms "global" and "local" are used, for which see the * accompanying paper. A brief summary of the 6 routines follows: * * BBINIT - initializes the system. * BBQR - the primary decomposition routine. It is called at least * once for each local set and generates the R matrix in * the QR decomposition of the design matrix. * BBSLVG - uses the R matrix to solve for the global solution vector. * BBSLVL - uses the R matrix to solve for one local vector per call. * BBCOVG - uses R to compute the global covariance matrix. * BBCOVL - uses R to compute the local and local-global covariances. * * These routines are generally called in the order listed. * * The routines make heavy use of the BLAS packages and require the * following routines: * * Level 1 BLAS: SCOPY, SASUM, SNRM2, SDOT, SAXPY * Level 2 BLAS: SGEMV, SGER, STRMV * Level 3 BLAS: SGEMM, SSYMM, STRMM, STRSM * LAPACK: SLARFG, STRTRI, SLAPY2, XERBLA * * LAPACK (indirect calls): ILAENV, SLAMCH, STRTI2, LSAME * * Written by: R. Ray Sept. 1992 * * SUBROUTINE BBINIT( NG, NRHS, RG, LRG, RMS ) * * Function - initializes system by zeroing global part of R matrix. * * Arguments - * * NG (input) INTEGER * Number of global parameters, and the order of the RG matrix. * (If NG is not known exactly at the beginning of processing, * then this number should be the largest number expected.) * * NRHS (input) INTEGER * The number of right-hand-side vectors in this problem. * * RG (output) REAL array, dimension (LRG,NG) * The global part of the R matrix. It is triangular and its * strictly lower triangle is never accessed. * * LRG (input) INTEGER * First dimension of RG in the calling program. * Must be >= NG. * * RMS (output) REAL array, dimension NRHS * Will contain the least-squares rms residuals, is initialized * to zero by this routine. * INTEGER NG, NRHS, LRG REAL RG(LRG,*), RMS(*), ZERO PARAMETER (ZERO=0.0) INTEGER I, J * * check size of LRG * ----------------- IF (NG.GT.LRG) THEN CALL XERBLA( 'BBINIT', 4 ) RETURN ENDIF * zero out RG and RMS * ------------------- DO 100 J=1,NG DO 100 I=1,J RG(I,J) = ZERO 100 CONTINUE DO 200 J=1,NRHS RMS(J) = ZERO 200 CONTINUE RETURN END *======================================================================= SUBROUTINE BBQR( NG,NL,NRHS, ! <-- problem size * NROWS,AG,LAG,AL,LAL,Y,LY, ! <-- design matrix * RL,LRL,RG,LRG,T,LT,CL,LCL,CG,LCG, ! <-- R matrix * WORK,RMS,NEWSET ) ! <-- miscellaneous * * Function - process NROWS of the design matrix to compute or update * the QR decomposition for a particular local set. * Must be called at least once for each local set. * * Arguments - * * NG (input) INTEGER * Number of global parameters; the order of the RG matrix. * * NL (input) INTEGER * Number of local parameters in the current set; the order * of the RL matrix. * * NRHS (input) INTEGER * The number of right-hand-side vectors in this problem. * * NROWS (input) INTEGER * Number of rows of AG and AL to process during this call. * * AG (input/output) REAL array, dimension (LAG,NROWS) * Rectangular submatrix of global parameter partials, stored * by rows. On exit, the data are destroyed. * * LAG (input) INTEGER * First dimension of AG in calling program. Must be >= NG. * * AL (input/output) REAL array, dimension (LAL,NROWS) * Submatrix of local parameter partials, stored by rows, for * use in computing R. On exit, the data are destroyed. * * LAL (input) INTEGER * First dimension of AL in calling program. Must be >= NL. * * Y (input/output) REAL array, dimension (LY,NRHS) * Right-hand-side vectors containing the observations for use * in computing CG and CL. On exit, the data are destroyed. * * LY (input) INTEGER * First dimension of Y in calling program. Must be >= NROWS. * * RL (input/output) REAL array, dimension (LRL,NL) * Triangular part of the local R matrix for the current set of * local parameters. Its strictly lower triangle is not used. * It is initialized to zero if NEWSET = TRUE. * * LRL (input) INTEGER * First dimension of RL in calling program. Must be >= NL. * * RG (input/output) REAL array, dimension (LRG,NG) * The global part of the R matrix. It is upper triangular; * its strictly lower triangle is not accessed. * * LRG (input) INTEGER * First dimension of RG in calling program. Must be >= NG. * * T (input/output) REAL array, dimension (LT,NG). * NL x NG array containing the border block of R for the * current local set. Is initialized to zero if NEWSET = TRUE. * * LT (input) INTEGER * First dimension of T in calling program. Must be >= NL. * * CL (input/output) REAL array, dimension (NCL,NRHS) * Right-hand-side vectors of decomposed observation vectors, * for the current set of local parameters. * * LCL (input) INTEGER * First dimension of CL in calling program. Must be >= NL. * * CG (input/output) REAL array, dimension (NCG,NRHS) * Right-hand-side vectors of decomposed global vectors. * * LCG (input) INTEGER * First dimension of CG in calling program. Must be >= NG. * * WORK (workspace) REAL array, minimum size 2*(NL + NG). * * RMS (input/output) REAL array, dimension NRHS * Contains the least-squares rms residuals according to * the observation equations already processed. * * NEWSET (input) LOGICAL * Logical flag set to .FALSE. if this call is to process * observations from a local set previously used, in which * case, RL and T are updated; otherwise RL & T are initialized. * I.e., NEWSET = .TRUE. if the rows of AL are coming from a * diagonal block not previously accessed. * * Usage note - * Usually, the arrays RG,RL,T,CG,CL are not modified by the driver, * but RL,T,CL must be saved for each local set if routines BBSLVL or * BBCOVL are to be called later. * INTEGER NG, NL, NRHS, NROWS, LAG, LAL, LY INTEGER LRL, LRG, LT, LCL, LCG REAL AG(LAG,*), AL(LAL,*), Y(LY,*) REAL RL(LRL,*), RG(LRG,*), T(LT,*), CL(LCL,*), CG(LCG,*) REAL WORK(*), RMS(*) LOGICAL NEWSET * REAL ZERO, ONE, TAU, Q11, TEMP, FACT INTEGER I, J, NN, JR, NR, IT, IR, IW PARAMETER (ZERO=0.0, ONE=1.0) REAL SASUM, SDOT, SNRM2, SLAPY2 EXTERNAL SASUM, SDOT, SNRM2, SLAPY2 * * If these data are for a new set of local parms, zero RL & T * ----------------------------------------------------------- IF (NEWSET) THEN DO 30 J=1,NL DO 10 I=1,J RL(I,J) = ZERO 10 CONTINUE DO 20 I=1,NG T(J,I) = ZERO 20 CONTINUE 30 CONTINUE ENDIF * * Eliminate AL one column at a time * --------------------------------- JR = 1 NR = NROWS DO 50 J=1,NL NN = NL - J + 1 IF (RL(J,J).EQ.ZERO) THEN * * If Jth row of RL is empty, fill it * ---------------------------------- IF ( SASUM(NN,RL(J,J),LRL).EQ.ZERO .AND. * SASUM(NG,T(J,1),LT).EQ.ZERO ) THEN CALL SCOPY( NN, AL(J,JR), 1, RL(J,J), LRL ) CALL SCOPY( NG, AG(1,JR), 1, T(J,1), LT ) CALL SCOPY( NRHS, Y(JR,1), LY, CL(J,1), LCL ) NR = NR - 1 JR = JR + 1 IF (NR.EQ.0) RETURN ENDIF ENDIF * * Determine Householder vector to eliminate Jth column of AL * ---------------------------------------------------------- CALL SLARFG( NR+1, RL(J,J), AL(J,JR), LAL, TAU ) IF (TAU.EQ.ZERO) GO TO 50 Q11 = ONE - TAU * Update RL, T, CL, & remainder of AL * ----------------------------------- IT = 1 IR = IT + NG IW = IR + NN - 1 CALL SCOPY( NG, T(J,1), LT, WORK(IT), 1 ) IF (J.LT.NL) THEN CALL SCOPY( NN-1, RL(J,J+1), LRL, WORK(IR), 1 ) * * Update Jth row of RL * -------------------- CALL SGEMV( 'N', NN-1, NR, -TAU, AL(J+1,JR), LAL, * AL(J,JR), LAL, Q11, RL(J,J+1), LRL ) * * Update remaining columns of AL * ------------------------------ CALL SGEMV( 'N', NN-1, NR, ONE, AL(J+1,JR), LAL, * AL(J,JR), LAL, ZERO, WORK(IW), 1 ) CALL SGER( NN-1, NR, -TAU, WORK(IW), 1, * AL(J,JR), LAL, AL(J+1,JR), LAL ) CALL SGER( NN-1, NR, -TAU, WORK(IR), 1, * AL(J,JR), LAL, AL(J+1,JR), LAL ) ENDIF * * Update Jth row of T * ------------------- CALL SGEMV( 'N', NG, NR, -TAU, AG(1,JR), LAG, * AL(J,JR), LAL, Q11, T(J,1), LT ) * * Update AG * --------- CALL SGEMV( 'N', NG, NR, ONE, AG(1,JR), LAG, * AL(J,JR), LAL, ZERO, WORK(IW), 1 ) CALL SGER( NG, NR, -TAU, WORK(IW), 1, * AL(J,JR), LAL, AG(1,JR), LAG ) CALL SGER( NG, NR, -TAU, WORK(IT), 1, * AL(J,JR), LAL, AG(1,JR), LAG ) * * Update right-hand-side vectors * ------------------------------ DO 40 I=1,NRHS TEMP = SDOT( NR, AL(J,JR), LAL, Y(JR,I), 1 ) FACT = -TAU*(CL(J,I) + TEMP) CL(J,I) = Q11*CL(J,I) - TAU*TEMP CALL SAXPY( NR, FACT, AL(J,JR), LAL, Y(JR,I), 1 ) 40 CONTINUE 50 CONTINUE * * Eliminate AG one column at a time * --------------------------------- DO 70 J=1,NG NN = NG - J + 1 IF (RG(J,J).EQ.ZERO) THEN * * If Jth row of RG is empty, fill it * ---------------------------------- IF ( SASUM(NN,RG(J,J),LRG).EQ.ZERO ) THEN CALL SCOPY( NN, AG(J,JR), 1, RG(J,J), LRG ) CALL SCOPY( NRHS, Y(JR,1), LY, CG(J,1), LCG ) NR = NR - 1 JR = JR + 1 IF (NR.EQ.0) RETURN ENDIF ENDIF * * Determine Householder vector to eliminate Jth column of AG * ---------------------------------------------------------- CALL SLARFG( NR+1, RG(J,J), AG(J,JR), LAG, TAU ) IF (TAU.EQ.ZERO) GO TO 70 Q11 = ONE - TAU * * Update RG & remainder of AL * --------------------------- IF (J.LT.NG) THEN IR = 1 IW = IR + NN - 1 CALL SCOPY( NN-1, RG(J,J+1), LRG, WORK(IR), 1 ) * * Update Jth row of RG * -------------------- CALL SGEMV( 'N', NN-1, NR, -TAU, AG(J+1,JR), LAG, * AG(J,JR), LAG, Q11, RG(J,J+1), LRG ) * * Update remaining columns of AG * ------------------------------ CALL SGEMV( 'N', NN-1, NR, ONE, AG(J+1,JR), LAG, * AG(J,JR), LAG, ZERO, WORK(IW), 1 ) CALL SGER( NN-1, NR, -TAU, WORK(IW), 1, * AG(J,JR), LAG, AG(J+1,JR), LAG ) CALL SGER( NN-1, NR, -TAU, WORK(IR), 1, * AG(J,JR), LAG, AG(J+1,JR), LAG ) ENDIF * * Update right-hand-side vectors * ------------------------------ DO 60 I=1,NRHS TEMP = SDOT( NR, AG(J,JR), LAG, Y(JR,I), 1 ) FACT = -TAU*(CG(J,I) + TEMP) CG(J,I) = Q11*CG(J,I) - TAU*TEMP CALL SAXPY( NR, FACT, AG(J,JR), LAG, Y(JR,I), 1 ) 60 CONTINUE 70 CONTINUE * * Update rms of residuals * ----------------------- DO 80 I=1,NRHS TEMP = SNRM2( NR, Y(JR,I), 1 ) RMS(I) = SLAPY2( RMS(I), TEMP ) 80 CONTINUE RETURN END *======================================================================= SUBROUTINE BBSLVG( NG, NRHS, RG, LRG, CG, LCG, INFO ) * * Function - solves for the global parameters after decomposition * has been completed. * * Arguments - * * NG (input) INTEGER * Number of global parameters, and the order of the RG matrix. * * NRHS (input) INTEGER * The number of right-hand-side vectors in this problem. * * RG (input) REAL array, dimension (LRG,NG) * The global part of the triangular R matrix as output by BBQR. * * LRG (input) INTEGER * First dimension of RG in the calling program. * * CG (input/output) REAL array, dimension (LCG,NRHS) * On entry, contains the right-hand-side vectors of the * global part of the decomposed design matrix. * On exit, contains the NRHS solution vectors for the * global parameters. * * LCG (input) INTEGER * First dimension of CG in calling program. Must be >= NG. * * INFO (output) INTEGER * Error flag; 0 denotes successful exit. * If > 0, a zero was detected along the diagonal of RG. * The matrix is singular and no solution is computed. * If < 0, an error was detected in the calling arguments. * * INTEGER LRG, LCG, NG, NRHS, INFO REAL RG(LRG,*), CG(LCG,*) REAL ZERO, ONE PARAMETER (ZERO=0.0, ONE=1.0) INTEGER I * * Check some of the arguments * --------------------------- INFO = 0 IF (NG.GT.LRG) THEN INFO = -4 ELSE IF (NG.GT.LCG) THEN INFO = -6 ENDIF IF (INFO.NE.0) THEN CALL XERBLA( 'BBSLVG', -INFO ) RETURN ENDIF * * Check that RG matrix is not singular * ------------------------------------ DO 10 I=1,NG IF (RG(I,I).EQ.ZERO) THEN INFO = I RETURN ENDIF 10 CONTINUE * * Compute solutions via level-3 BLAS routine * ------------------------------------------ IF (NRHS.LT.1) RETURN CALL STRSM( 'L', 'U', 'N', 'N', NG, NRHS, ONE, RG, LRG, * CG, LCG ) * RETURN END *======================================================================= SUBROUTINE BBSLVL( NG ,NL, NRHS, RL, LRL, T, LT, * CL, LCL, XG, LXG, INFO ) * * Function - solves for the local parameters for one local set. * * Arguments - * * NG (input) INTEGER * Number of global parameters; the order of the RG matrix. * * NL (input) INTEGER * Number of local parameters in this set; the order of RL. * * NRHS (input) INTEGER * The number of right-hand-side vectors in this problem. * * RL (input) REAL array, dimension (LRL,NL) * The local part of the triangular R matrix for one set of * local parameters, as output by BBQR. * * LRL (input) INTEGER * First dimension of RL in calling program. Must be >= NL. * * T (input) REAL array, dimension (LT,NG) * Contains one border block in the R matrix, corresponding to * the current set of local parameters, as output by BBQR. * * LT (input) INTEGER * First dimension of T in calling program. Must be >= NL. * * CL (input/output) REAL array, dimension (LCL,NRHS) * On entry, contains the right-hand-side vectors of one * local part of the decomposed observation equations. * On exit, contains the NRHS solution vectors for the * current local parameters. * * LCL (input) INTEGER * First dimension of CL in calling program. Must be >= NL. * * XG (input) REAL array, dimension (LXG,NRHS) * Contains the global solution vectors, as returned by BBSLVG. * * LXG (input) INTEGER * First dimension of XG in calling program. Must be >= NG. * * INFO (output) INTEGER * Error flag; 0 denotes successful exit. * If > 0, a zero was detected along the diagonal of RL. * The matrix is singular and no solution was computed. * If < 0, an error was detected in the calling arguments. * * INTEGER NG, NL, NRHS, LRL, LT, LCL, LXG, INFO REAL RL(LRL,*), T(LT,*), CL(LCL,*), XG(LXG,*) REAL ZERO, ONE PARAMETER (ZERO=0.0, ONE=1.0) INTEGER I * * Check some of the arguments * --------------------------- INFO = 0 IF (NL.GT.LRL) THEN INFO = -5 ELSE IF (NL.GT.LT) THEN INFO = -7 ELSE IF (NL.GT.LCL) THEN INFO = -9 ELSE IF (NG.GT.LXG) THEN INFO = -11 ENDIF IF (INFO.NE.0) THEN CALL XERBLA( 'BBSLVG', -INFO ) RETURN ENDIF * * Check that RL matrix is not singular * ------------------------------------ IF (NL.EQ.0) RETURN DO 10 I=1,NL IF (RL(I,I).EQ.ZERO) THEN INFO = I RETURN ENDIF 10 CONTINUE * * Subtract T*XG from CL * --------------------- CALL SGEMM( 'N', 'N', NL, NRHS, NG, -ONE, T, LT, XG, LXG, * ONE, CL, LCL ) * * Solve RL*XL = CL by back substitution * ------------------------------------- CALL STRSM( 'L', 'U', 'N', 'N', NL, NRHS, ONE, RL, LRL, * CL, LCL ) * RETURN END *======================================================================= SUBROUTINE BBCOVG( NG, RG, LRG, WORK, INFO ) * * Function - computes the global covariance matrix. * * Arguments - * * NG (input) INTEGER * Number of global parameters; the order of the RG matrix. * * RG (input/output) REAL array, dimension (LRG,NG) * On entry, RG contains the global part of the R matrix. * On exit, it contains the upper triangle of the symmetric * covariance matrix for the global parameters. * The strictly lower triangle of RG is never accessed. * * LRG (input) INTEGER * First dimension of RG in the calling program. * * WORK (workspace) REAL array, minimum size NG. * * INFO (output) INTEGER * Error flag; 0 denotes successful exit. * If > 0, a zero was detected along the diagonal of RG. * The matrix is singular and no covariance is computed. * If < 0, an error was detected in the calling arguments. * INTEGER NG, LRG, INFO REAL RG(LRG,*), WORK(*) INTEGER I, NN * * Compute inverse of RG * --------------------- CALL STRTRI( 'U', 'N', NG, RG, LRG, INFO ) IF (INFO.NE.0) RETURN * * Compute inverse times its transpose * ----------------------------------- DO 100 I=1,NG NN = NG - I + 1 CALL SCOPY( NN, RG(I,I), LRG, WORK, 1 ) CALL STRMV( 'U', 'N', 'N', NN, RG(I,I), LRG, WORK, 1 ) CALL SCOPY( NN, WORK, 1, RG(I,I), LRG ) 100 CONTINUE * RETURN END *======================================================================= SUBROUTINE BBCOVL( NG, NL, RL, LRL, T, LT, CVG, LCVG, WORK, INFO ) * * Function - computes the local and local-global covariances for one * local set. * * Arguments - * * NG (input) INTEGER * Number of global parameters, and the order of the RG matrix. * * NL (input) INTEGER * Number of local parameters in this set; the order of RL. * * RL (input/output) REAL array, dimension (LRL,NL) * On entry, contains the upper triangular local-parameter * part of the R matrix, as output by BBQR, for this set of * local parameters. * On exit, it contains the upper triangle of the symmetric * covariance matrix for this same set of local parameters. * The strictly lower triangle of RL is never accessed. * * LRL (input) INTEGER * First dimension of RL in calling program. Must be >= NL. * * T (input/output) REAL array, dimension (LT,NG) * On entry, contains one border block in the R matrix, * corresponding to this set of local parameters, as output * by BBQR. * On exit, contains the local-global covariances for this * local set. * * LT (input) INTEGER * First dimension of T in calling program. Must be >= NL. * * CVG (input) REAL array, dimension (LCVG,NG) * Upper triangular matrix of global covariances, as output * by BBCOVG. * * LCVG (input) INTEGER * First dimension of CVG in calling program. * Must be >= NG. * * WORK (workspace) REAL array, minimum size (NG + 1)*NL. * * INFO (output) INTEGER * Error flag; 0 denotes successful exit. * If > 0, a zero was detected along the diagonal of RL. * The matrix is singular and no covariances are computed. * If < 0, an error was detected in the calling arguments. * * INTEGER NG, NL, LRL, LT, LCVG, INFO REAL RL(LRL,*), T(LT,*), CVG(LCVG,*), WORK(*) REAL ZERO, ONE PARAMETER (ZERO=0.0, ONE=1.0) INTEGER I, J, NN, IW * * Compute inverse of RL * --------------------- CALL STRTRI( 'U', 'N', NL, RL, LRL, INFO ) IF (INFO.NE.0) RETURN * * Compute Inv(RL) * T and store in WORK * ------------------------------------- DO 10 J=1,NG CALL SCOPY( NL, T(1,J), 1, WORK((J-1)*NL+1), 1 ) 10 CONTINUE CALL STRMM( 'L', 'U', 'N', 'N', NL, NG, ONE, RL, LRL, WORK, NL ) * * Compute global-local cross-covariances and store in T * ----------------------------------------------------- CALL SSYMM( 'R', 'U', NL, NG, -ONE, CVG, LCVG, WORK, NL, * ZERO, T, LT ) * * Compute local covariance matrix row by row * ------------------------------------------ IW = NG*NL + 1 DO 20 I=1,NL NN = NL - I + 1 CALL SCOPY( NN, RL(I,I), LRL, WORK(IW), 1 ) CALL STRMV( 'U', 'N', 'N', NN, RL(I,I), LRL, WORK(IW), 1 ) CALL SCOPY( NN, WORK(IW), 1, RL(I,I), LRL ) CALL SGEMV( 'N', NN, NG, -ONE, WORK(I), NL, T(I,1), LT, * ONE, RL(I,I), LRL ) 20 CONTINUE * RETURN END c*** bbtest.f PROGRAM BBTEST * * Example test driver for BBQR * * This program solves a system of 20 equations (via least squares) * for 3 global parameters and 2 sets of local parameters, each set * containing 2 parameters. * There are 10 observation equations corresponding to each local set. * There are 2 right-hand-side vectors, one being an exact solution. * First the system is solved using straightforward calls to * LAPACK routines for a full (dense) matrix. These results * may then be compared to the BBQR,BBSLVx,BBCOVx results. * * R. Ray Oct 1992 * * NP - number of observation equations * N - number of unknowns for full-matrix case * NG - number of global parameters to solve for * NL - maximum number of local parameters in any one set * ------------------------------------------------------ PARAMETER (NP=20,N=7,NG=3,NL=2,NRHS=2) * arrays for full-matrix computations: * ------------------------------------ PARAMETER (LWORK=NP*64+2*NP+2*NG) DIMENSION A(NP,N),Y(NP,NRHS),X(N),P(N,NP),WORK(LWORK),B(NP,NRHS) DIMENSION COV(N,N) * arrays for BBQR global data * --------------------------- DIMENSION CG(NG,NRHS),AG(NG,NP),RG(NG,NG) * arrays for BBQR local data * (last index is for 2 local sets) * -------------------------------- DIMENSION AL(NL,NP/2,2),CL(NL,NRHS,2),RL(NL,NL,2) DIMENSION TGL(NL,NG,2),RMS(2) * * True solution: DATA X/ 1.0, -0.5, 0.5, -2.0, -1.2, 1.2, 1.0/ * * Design matrix of partial derivatives * (10 equations for each local set, with total of 7 unknowns) * ------------------------------------ DATA P/ 1., 1., 0., 0., 0., 0., 0., * 0., 1., 0., 0., 1., 0., 0., * 2., -1., 0., 0., 1., 1.,-1., * 2., 2., 0., 0., 2., 1., 1., * 3., 2., 0., 0., 2.,-1., 0., * -1., -1., 0., 0., 2.,-2., 1., * -1., 0., 0., 0., -2.,-2., 1., * -2., -2., 0., 0., -2.,-3., 3., * -2., -1., 0., 0., 3., 5., 2., * 1., -1., 0., 0., 4., 1., 3., * 0., 0., 0., 0., 0., 1., 0., * 0., 0., 0., 1., 1., 0., 0., * 0., 0., 2., 0., 1., 1.,-1., * 0., 0., 2.,-1., 2., 1., 1., * 0., 0., 1., 1., 2.,-1., 0., * 0., 0., -1., 3., 2.,-2., 1., * 0., 0., 3., 0., -2.,-2., 1., * 0., 0., -3.,-1., -2.,-3., 3., * 0., 0., -2., 2., 3., 5., 2., * 0., 0., -2.,-2., 4., 1., 3. / * Build & solve full-matrix system * -------------------------------- DO 40 J=1,NP B(J,1) = SDOT(N,P(1,J),1,X,1) B(J,2) = B(J,1) + 0.1 DO 30 I=1,N A(J,I) = P(I,J) 30 CONTINUE CALL SSYR( 'U',N,1.,P(1,J),1,COV,N ) 40 CONTINUE CALL SGELS( 'N',NP,N,NRHS,A,NP,B,NP,WORK,LWORK,INFO ) RMS(1) = SNRM2( NP-N, B(N+1,1), 1 ) RMS(2) = SNRM2( NP-N, B(N+1,2), 1 ) CALL SPOTRF( 'U',N,COV,N,INFO ) CALL SPOTRI( 'U',N,COV,N,INFO ) WRITE(6,50) 50 FORMAT(/' True R'' matrix:'/) DO 60 J=1,N 60 WRITE(6,100) (A(I,J),I=1,J) WRITE(6,65) 65 FORMAT(/' True covariance:'/) DO 66 J=1,N 66 WRITE(6,100) (COV(I,J),I=1,J) WRITE(6,70) WRITE(6,75) (B(I,1),I=1,N) WRITE(6,75) (B(I,2),I=1,N) 70 FORMAT(/' True solution vectors:') 75 FORMAT(/1X,8F10.3) WRITE(6,80) RMS 80 FORMAT(/' True residual rms:',F10.3,' and',F10.3) 100 FORMAT(1X,8F10.3) * Build BBQR 'observations' * ------------------------- DO 120 J=1,NP Y(J,1) = SDOT(N,P(1,J),1,X,1) Y(J,2) = Y(J,1) + 0.1 DO 110 I=1,NG AG(I,J) = P(2*NL+I,J) 110 CONTINUE 120 CONTINUE DO 140 J=1,NP/2 DO 130 I=1,NL AL(I,J,1) = P(I,J) AL(I,J,2) = P(I+NL,J+NP/2) 130 CONTINUE 140 CONTINUE * Test BBQR on same data * ----------------------- *-----Initialize the system CALL BBINIT( NG, NRHS, RG, NG, RMS ) *-----Process all 10 equations for 1st local set CALL BBQR( NG,NL,NRHS, 10,AG,NG,AL,NL,Y,NP, * RL,NL,RG,NG,TGL,NL,CL,NL,CG,NG, * WORK,RMS,.TRUE. ) WRITE(6,*) 'Set 1 processed.' *-----Process 10 equations for 2nd local set, one equation at a time CALL BBQR( NG,NL,NRHS, 1,AG(1,11),NG,AL(1,1,2),NL,Y(11,1),NP, * RL(1,1,2),NL,RG,NG,TGL(1,1,2),NL,CL(1,1,2),NL,CG,NG, * WORK,RMS,.TRUE. ) DO 200 J=2,10 CALL BBQR(NG,NL,NRHS,1,AG(1,10+J),NG,AL(1,J,2),NL,Y(10+J,1),NP, * RL(1,1,2),NL,RG,NG,TGL(1,1,2),NL,CL(1,1,2),NL,CG,NG, * WORK,RMS,.FALSE. ) 200 CONTINUE WRITE(6,*) 'Set 2 processed.' WRITE(6,210) 210 FORMAT(/' Results for BBQR:'//' R'' matrix:') WRITE(6,*) 'Global:' DO 220 J=1,NG 220 WRITE(6,100) (RG(I,J),I=1,J) DO 230 K=1,2 WRITE(6,222) K 222 FORMAT(' Local set',I3,':') DO 225 J=1,NL 225 WRITE(6,100) (RL(I,J,K),I=1,J) DO 227 J=1,NG 227 WRITE(6,100) (TGL(I,J,K),I=1,NL) 230 CONTINUE * Compute global & local solution vectors * --------------------------------------- CALL BBSLVG( NG, NRHS, RG, NG, CG, NG, INFO ) CALL BBSLVL( NG,NL,NRHS,RL(1,1,1),NL,TGL(1,1,1),NL, * CL(1,1,1),NL,CG,NG,INFO ) CALL BBSLVL( NG,NL,NRHS,RL(1,1,2),NL,TGL(1,1,2),NL, * CL(1,1,2),NL,CG,NG,INFO ) WRITE(6,240) WRITE(6,75) (CG(I,1),I=1,NG) WRITE(6,75) (CG(I,2),I=1,NG) 240 FORMAT(/' Global solution vectors:') WRITE(6,270) WRITE(6,75) CL(1,1,1),CL(2,1,1) WRITE(6,75) CL(1,2,1),CL(2,2,1) 270 FORMAT(/' Solution vectors for local set #1:') WRITE(6,290) WRITE(6,75) CL(1,1,2),CL(2,1,2) WRITE(6,75) CL(1,2,2),CL(2,2,2) 290 FORMAT(/' Solution vectors for local set #2:') WRITE(6,245) RMS 245 FORMAT(/' Residual rms:',2F10.3) * Compute covariances * ------------------- CALL BBCOVG( NG, RG, NG, WORK, INFO ) CALL BBCOVL( NG, NL, RL, NL, TGL, NL, RG, NG, WORK, INFO ) CALL BBCOVL( NG, NL, RL(1,1,2), NL, TGL(1,1,2), NL, RG, NG, * WORK, INFO ) WRITE(6,250) 250 FORMAT(/' Global covariance matrix:'/) DO 260 J=1,NG 260 WRITE(6,100) (RG(I,J),I=1,J) DO 300 K=1,2 WRITE(6,280) K 280 FORMAT(/' Local covariance matrices for set #',I1,':'/) WRITE(6,100) RL(1,1,K) WRITE(6,100) RL(1,2,K),RL(2,2,K) DO 293 J=1,NG 293 WRITE(6,100) TGL(1,J,K),TGL(2,J,K) 300 CONTINUE STOP END c*** support.f real function sasum(n,sx,incx) c c takes the sum of the absolute values. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c real sx(1),stemp integer i,incx,m,mp1,n,nincx c sasum = 0.0e0 stemp = 0.0e0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx stemp = stemp + abs(sx(i)) 10 continue sasum = stemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + abs(sx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) * + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5)) 50 continue 60 sasum = stemp return end subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c real sx(1),sy(1),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0.0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sy(iy) + sa*sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sy(i) + sa*sx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = sy(i) + sa*sx(i) sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end subroutine scopy(n,sx,incx,sy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c real sx(1),sy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 sy(i) = sx(i) sy(i + 1) = sx(i + 1) sy(i + 2) = sx(i + 2) sy(i + 3) = sx(i + 3) sy(i + 4) = sx(i + 4) sy(i + 5) = sx(i + 5) sy(i + 6) = sx(i + 6) 50 continue return end real function sdot(n,sx,incx,sy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c real sx(1),sy(1),stemp integer i,incx,incy,ix,iy,m,mp1,n c stemp = 0.0e0 sdot = 0.0e0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + sx(i)*sy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end real function snrm2 ( n, sx, incx) integer next real sx(1), cutlo, cuthi, hitest, sum, xmax, zero, one data zero, one /0.0e0, 1.0e0/ c c euclidean norm of the n-vector stored in sx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of sqrt(u/eps) over all known machines. c cuthi = minimum of sqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 4.441e-16, 1.304e19 / c if(n .gt. 0) go to 10 snrm2 = zero go to 300 c 10 assign 30 to next sum = zero nn = n * incx c begin main loop i = 1 20 go to next,(30, 50, 70, 110) 30 if( abs(sx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( sx(i) .eq. zero) go to 200 if( abs(sx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 i = j assign 110 to next sum = (sum / sx(i)) / sx(i) 105 xmax = abs(sx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( abs(sx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( abs(sx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / sx(i))**2 xmax = abs(sx(i)) go to 200 c 115 sum = sum + (sx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j =i,nn,incx if(abs(sx(j)) .ge. hitest) go to 100 95 sum = sum + sx(j)**2 snrm2 = sqrt( sum ) go to 300 c 200 continue i = i + incx if ( i .le. nn ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c snrm2 = xmax * sqrt(sum) 300 continue return end subroutine sscal(n,sa,sx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to 1. c jack dongarra, linpack, 3/11/78. c real sa,sx(1) integer i,incx,m,mp1,n,nincx c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx sx(i) = sa*sx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m sx(i) = sa*sx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 sx(i) = sa*sx(i) sx(i + 1) = sa*sx(i + 1) sx(i + 2) = sa*sx(i + 2) sx(i + 3) = sa*sx(i + 3) sx(i + 4) = sa*sx(i + 4) 50 continue return end SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. REAL ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. REAL A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * SGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - REAL array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. Local Scalars .. REAL TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of SGEMV . * END SUBROUTINE SGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. REAL ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. REAL A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * SGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. Local Scalars .. REAL TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of SGER . * END SUBROUTINE SSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * .. Scalar Arguments .. REAL ALPHA INTEGER INCX, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. REAL A( LDA, * ), X( * ) * .. * * Purpose * ======= * * SSYR performs the symmetric rank 1 operation * * A := alpha*x*x' + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of A is not referenced. On exit, the * upper triangular part of the array A is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of A is not referenced. On exit, the * lower triangular part of the array A is overwritten by the * lower triangular part of the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. Local Scalars .. REAL TEMP INTEGER I, INFO, IX, J, JX, KX * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SSYR ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Set the start point in X if the increment is not unity. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * IF( LSAME( UPLO, 'U' ) )THEN * * Form A when A is stored in upper triangle. * IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) DO 10, I = 1, J A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = KX DO 30, I = 1, J A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JX = JX + INCX 40 CONTINUE END IF ELSE * * Form A when A is stored in lower triangle. * IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) DO 50, I = J, N A( I, J ) = A( I, J ) + X( I )*TEMP 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = JX DO 70, I = J, N A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF * RETURN * * End of SSYR . * END SUBROUTINE STRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. REAL A( LDA, * ), X( * ) * .. * * Purpose * ======= * * STRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. Local Scalars .. REAL TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of STRMV . * END SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC REAL ALPHA, BETA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * SGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - REAL array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of SGEMM . * END SUBROUTINE SSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO INTEGER M, N, LDA, LDB, LDC REAL ALPHA, BETA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * SSYMM performs one of the matrix-matrix operations * * C := alpha*A*B + beta*C, * * or * * C := alpha*B*A + beta*C, * * where alpha and beta are scalars, A is a symmetric matrix and B and * C are m by n matrices. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether the symmetric matrix A * appears on the left or right in the operation as follows: * * SIDE = 'L' or 'l' C := alpha*A*B + beta*C, * * SIDE = 'R' or 'r' C := alpha*B*A + beta*C, * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the symmetric matrix A is to be * referenced as follows: * * UPLO = 'U' or 'u' Only the upper triangular part of the * symmetric matrix is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of the * symmetric matrix is to be referenced. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix C. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix C. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, ka ), where ka is * m when SIDE = 'L' or 'l' and is n otherwise. * Before entry with SIDE = 'L' or 'l', the m by m part of * the array A must contain the symmetric matrix, such that * when UPLO = 'U' or 'u', the leading m by m upper triangular * part of the array A must contain the upper triangular part * of the symmetric matrix and the strictly lower triangular * part of A is not referenced, and when UPLO = 'L' or 'l', * the leading m by m lower triangular part of the array A * must contain the lower triangular part of the symmetric * matrix and the strictly upper triangular part of A is not * referenced. * Before entry with SIDE = 'R' or 'r', the n by n part of * the array A must contain the symmetric matrix, such that * when UPLO = 'U' or 'u', the leading n by n upper triangular * part of the array A must contain the upper triangular part * of the symmetric matrix and the strictly lower triangular * part of A is not referenced, and when UPLO = 'L' or 'l', * the leading n by n lower triangular part of the array A * must contain the lower triangular part of the symmetric * matrix and the strictly upper triangular part of A is not * referenced. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, n ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - REAL array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n updated * matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, K, NROWA REAL TEMP1, TEMP2 * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Set NROWA as the number of rows of A. * IF( LSAME( SIDE, 'L' ) )THEN NROWA = M ELSE NROWA = N END IF UPPER = LSAME( UPLO, 'U' ) * * Test the input parameters. * INFO = 0 IF( ( .NOT.LSAME( SIDE, 'L' ) ).AND. $ ( .NOT.LSAME( SIDE, 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO, 'L' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 9 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 12 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SSYMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( LSAME( SIDE, 'L' ) )THEN * * Form C := alpha*A*B + beta*C. * IF( UPPER )THEN DO 70, J = 1, N DO 60, I = 1, M TEMP1 = ALPHA*B( I, J ) TEMP2 = ZERO DO 50, K = 1, I - 1 C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) TEMP2 = TEMP2 + B( K, J )*A( K, I ) 50 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ TEMP1*A( I, I ) + ALPHA*TEMP2 END IF 60 CONTINUE 70 CONTINUE ELSE DO 100, J = 1, N DO 90, I = M, 1, -1 TEMP1 = ALPHA*B( I, J ) TEMP2 = ZERO DO 80, K = I + 1, M C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) TEMP2 = TEMP2 + B( K, J )*A( K, I ) 80 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ TEMP1*A( I, I ) + ALPHA*TEMP2 END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form C := alpha*B*A + beta*C. * DO 170, J = 1, N TEMP1 = ALPHA*A( J, J ) IF( BETA.EQ.ZERO )THEN DO 110, I = 1, M C( I, J ) = TEMP1*B( I, J ) 110 CONTINUE ELSE DO 120, I = 1, M C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J ) 120 CONTINUE END IF DO 140, K = 1, J - 1 IF( UPPER )THEN TEMP1 = ALPHA*A( K, J ) ELSE TEMP1 = ALPHA*A( J, K ) END IF DO 130, I = 1, M C( I, J ) = C( I, J ) + TEMP1*B( I, K ) 130 CONTINUE 140 CONTINUE DO 160, K = J + 1, N IF( UPPER )THEN TEMP1 = ALPHA*A( J, K ) ELSE TEMP1 = ALPHA*A( K, J ) END IF DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP1*B( I, K ) 150 CONTINUE 160 CONTINUE 170 CONTINUE END IF * RETURN * * End of SSYMM . * END SUBROUTINE STRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL ALPHA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * STRMM performs one of the matrix-matrix operations * * B := alpha*op( A )*B, or B := alpha*B*op( A ), * * where alpha is a scalar, B is an m by n matrix, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) multiplies B from * the left or right as follows: * * SIDE = 'L' or 'l' B := alpha*op( A )*B. * * SIDE = 'R' or 'r' B := alpha*B*op( A ). * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B, and on exit is overwritten by the * transformed matrix. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE ELSE DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*A. * IF( UPPER )THEN DO 180, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE ELSE DO 220, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE END IF 210 CONTINUE 220 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 260, K = 1, N DO 240, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE END IF 240 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE ELSE DO 300, K = N, 1, -1 DO 280, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF 300 CONTINUE END IF END IF END IF * RETURN * * End of STRMM . * END SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL ALPHA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * STRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of STRSM . * END SUBROUTINE SSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC REAL ALPHA, BETA * .. Array Arguments .. REAL A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * SSYRK performs one of the symmetric rank k operations * * C := alpha*A*A' + beta*C, * * or * * C := alpha*A'*A + beta*C, * * where alpha and beta are scalars, C is an n by n symmetric matrix * and A is an n by k matrix in the first case and a k by n matrix * in the second case. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array C is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of C * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of C * is to be referenced. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. * * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. * * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with TRANS = 'N' or 'n', K specifies the number * of columns of the matrix A, and on entry with * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number * of rows of the matrix A. K must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - REAL array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN INFO = 2 ELSE IF( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SSYRK ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) )THEN * * Form C := alpha*A*A' + beta*C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of SSYRK . * END SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * SLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function SLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = SLAMC3( A, ONE ) C = SLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = SLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = SLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = SLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = SLAMC3( B / 2, -B / 100 ) C = SLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = SLAMC3( B / 2, B / 100 ) C = SLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = SLAMC3( B / 2, A ) T2 = SLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = SLAMC3( A, ONE ) C = SLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of SLAMC1 * END * ************************************************************************ * SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T REAL EPS, RMAX, RMIN * .. * * Purpose * ======= * * SLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) REAL * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) REAL * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) REAL * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. External Subroutines .. EXTERNAL SLAMC1, SLAMC4, SLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function SLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL SLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = SLAMC3( B, -HALF ) THIRD = SLAMC3( SIXTH, SIXTH ) B = SLAMC3( THIRD, -HALF ) B = SLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = SLAMC3( HALF, -C ) B = SLAMC3( HALF, C ) C = SLAMC3( HALF, -B ) B = SLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = SLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = SLAMC3( ONE, SMALL ) CALL SLAMC4( NGPMIN, ONE, LBETA ) CALL SLAMC4( NGNMIN, -ONE, LBETA ) CALL SLAMC4( GPMIN, A, LBETA ) CALL SLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine SLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute *