 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION DQDOTA(N,DB,QC,DX,INCX,DY,INCY) 00002 C D.P. DOT PRODUCT WITH EXTENDED PRECISION ACCUMULATION (AND RESULT) 00003 C QC AND DQDOTA ARE SET = DB + QC + SUM FOR I = 0 TO N-1 OF 00004 C DX(LX+I*INCX) * DY(LY+I*INCY), WHERE QC IS AN EXTENDED 00005 C PRECISION RESULT PREVIOUSLY COMPUTED BY DQDOTI OR DQDOTA 00006 C AND LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS 00007 C DEFINED IN A SIMILAR WAY USING INCY. THE MP PACKAGE BY 00008 C RICHARD P. BRENT IS USED FOR THE EXTENDED PRECISION ARITHMETIC. 00009 C 00010 C FRED T. KROGH, JPL, 1977, JUNE 1 00011 C2 00012 DOUBLE PRECISION DX(1), DY(1), DB 00013 INTEGER QC(10), QX(10), QY(10) 00014 C THE COMMON BLOCK FOR THE MP PACKAGE (MODIFIED TO GIVE IT A NAME) 00015 COMMON /MPCOM/ MPB, MPT, MPM, MPLUN, MPMXR, MPR(12) 00016 DATA I1 / 0 / 00017 C IF I1 IS 0 THE MP PACKAGE MUST BE INITIALIZED (MPBLAS SETS I1 = 1) 00018 IF (I1 .EQ. 0) CALL MPBLAS(I1) 00019 IF (DB .EQ. 0.D0) GO TO 20 00020 CALL MPCDM(DB, QX) 00021 CALL MPADD(QC, QX, QC) 00022 20 IF (N .EQ. 0) GO TO 40 00023 IX = 1 00024 IY = 1 00025 IF (INCX .LT. 0) IX = (-N + 1) * INCX + 1 00026 IF (INCY .LT. 0) IY = (-N + 1) * INCY + 1 00027 DO 30 I = 1,N 00028 CALL MPCDM(DX(IX), QX) 00029 CALL MPCDM(DY(IY), QY) 00030 CALL MPMUL(QX, QY, QX) 00031 CALL MPADD(QC, QX, QC) 00032 IX = IX + INCX 00033 IY = IY + INCY 00034 30 CONTINUE 00035 40 CALL MPCMD(QC, DQDOTA) 00036 RETURN 00037 END 00038 SUBROUTINE DTEST(LEN,DCOMP,DTRUE,DSIZE,DFAC) 00039 C1 ********************************* DTEST ************************** 00040 C 00041 C THIS SUBR COMPARES ARRAYS DCOMP() AND DTRUE() OF LENGTH LEN TO 00042 C SEE IF THE TERM BY TERM DIFFERENCES, MULTIPLIED BY DFAC, ARE 00043 C NEGLIGIBLE. 00044 C 00045 C C. L. LAWSON, JPL, 1974 DEC 10 00046 C2 00047 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00048 LOGICAL PASS 00049 DOUBLE PRECISION DCOMP(LEN),DTRUE(LEN),DSIZE(LEN),DFAC,DDIFF,DD 00050 C 00051 DO 10 I=1,LEN 00052 DD = DCOMP(I)-DTRUE(I) 00053 IF(DDIFF(DABS(DSIZE(I))+DABS(DFAC*DD),DABS(DSIZE(I))) .EQ. 0.D0) 00054 * GO TO 10 00055 C 00056 C HERE DCOMP(I) IS NOT CLOSE TO DTRUE(I). 00057 C 00058 IF(.NOT. PASS) GO TO 5 00059 C PRINT FAIL MESSAGE AND HEADER. 00060 PASS = .FALSE. 00061 WRITE(NPRINT,1000) 00062 WRITE(NPRINT,1001) 00063 5 WRITE(NPRINT,1002) ICASE,N,INCX,INCY,MODE,I, 00064 * DCOMP(I),DTRUE(I),DD,DSIZE(I) 00065 10 CONTINUE 00066 RETURN 00067 1000 FORMAT(1H+,39X,4HFAIL) 00068 1001 FORMAT(26H0CASE N INCX INCY MODE I, 00069 1 29X,7HCOMP(I),29X,7HTRUE(I),2X,10HDIFFERENCE, 00070 2 5X,7HSIZE(I)/1X) 00071 1002 FORMAT(1X,I4,I3,3I5,I3,2D36.18,2D12.4) 00072 END 00073 FUNCTION SDIFF(SA,SB) 00074 C1 ********************************* SDIFF ************************** 00075 C COMPUTES DIFFERENCE OF TWO NUMBERS. C. L. LAWSON, JPL 1974 FEB 15 00076 C2 00077 SDIFF=SA-SB 00078 RETURN 00079 END 00080 DOUBLE PRECISION FUNCTION DDIFF(DA,DB) 00081 C1 ********************************* DDIFF ************************** 00082 C COMPUTES DIFFERENCE OF TWO NUMBERS. C. L. LAWSON, JPL 1974 FEB 15 00083 C2 00084 DOUBLE PRECISION DA,DB 00085 DDIFF=DA-DB 00086 RETURN 00087 END 00088 SUBROUTINE STEST1(SCOMP1, STRUE1, SSIZE, SFAC) 00089 C1 ************************* STEST1 ***************************** 00090 C 00091 C THIS IS AN INTERFACE SUBROUTINE TO ACCOMODATE THE FORTRAN 00092 C REQUIREMENT THAT WHEN A DUMMY ARGUMENT IS AN ARRAY, THE 00093 C ACTUAL ARGUMENT MUST ALSO BE AN ARRAY OR AN ARRAY ELEMENT. 00094 C 00095 C C.L. LAWSON, JPL, 1978 DEC 6 00096 C2 00097 REAL SCOMP(1),STRUE(1),SSIZE(1) 00098 C 00099 SCOMP(1) = SCOMP1 00100 STRUE(1) = STRUE1 00101 CALL STEST(1,SCOMP,STRUE,SSIZE,SFAC) 00102 C 00103 RETURN 00104 END 00105 SUBROUTINE CTEST(LEN,CCOMP,CTRUE,CSIZE,SFAC) 00106 C1 **************************** CTEST ***************************** 00107 C 00108 C C.L. LAWSON, JPL, 1978 DEC 6 00109 C2 00110 COMPLEX CCOMP(LEN),CTRUE(LEN),CSIZE(LEN) 00111 REAL SFAC 00112 REAL SCOMP(20),STRUE(20),SSIZE(20) 00113 C 00114 DO 10 I=1,LEN 00115 SCOMP(2*I-1) = REAL(CCOMP(I)) 00116 SCOMP(2*I) = AIMAG(CCOMP(I)) 00117 STRUE(2*I-1) = REAL(CTRUE(I)) 00118 STRUE(2*I) = AIMAG(CTRUE(I)) 00119 SSIZE(2*I-1) = REAL(CSIZE(I)) 00120 SSIZE(2*I) = AIMAG(CSIZE(I)) 00121 10 CONTINUE 00122 C 00123 CALL STEST(2*LEN,SCOMP,STRUE,SSIZE,SFAC) 00124 RETURN 00125 END 00126 SUBROUTINE DTEST1(DCOMP1, DTRUE1, DSIZE, DFAC) 00127 C1 ************************* DTEST1 ***************************** 00128 C 00129 C THIS IS AN INTERFACE SUBROUTINE TO ACCOMODATE THE FORTRAN 00130 C REQUIREMENT THAT WHEN A DUMMY ARGUMENT IS AN ARRAY, THE 00131 C ACTUAL ARGUMENT MUST ALSO BE AN ARRAY OR AN ARRAY ELEMENT. 00132 C 00133 C C.L. LAWSON, JPL, 1978 DEC 6 00134 C2 00135 DOUBLE PRECISION DCOMP1, DTRUE1, DFAC 00136 DOUBLE PRECISION DCOMP(1),DTRUE(1),DSIZE(1) 00137 C 00138 DCOMP(1) = DCOMP1 00139 DTRUE(1) = DTRUE1 00140 CALL DTEST(1,DCOMP,DTRUE,DSIZE,DFAC) 00141 C 00142 RETURN 00143 END 00144 SUBROUTINE STEST(LEN,SCOMP,STRUE,SSIZE,SFAC) 00145 C1 ********************************* STEST ************************** 00146 C 00147 C THIS SUBR COMPARES ARRAYS SCOMP() AND STRUE() OF LENGTH LEN TO 00148 C SEE IF THE TERM BY TERM DIFFERENCES, MULTIPLIED BY SFAC, ARE 00149 C NEGLIGIBLE. 00150 C 00151 C C. L. LAWSON, JPL, 1974 DEC 10 00152 C2 00153 REAL SCOMP(LEN),STRUE(LEN),SSIZE(LEN),SFAC,SDIFF,SD 00154 LOGICAL PASS 00155 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00156 C 00157 DO 10 I=1,LEN 00158 SD = SCOMP(I)-STRUE(I) 00159 IF( SDIFF(ABS(SSIZE(I))+ABS(SFAC*SD), ABS(SSIZE(I))) .EQ. 0.) 00160 * GO TO 10 00161 C 00162 C HERE SCOMP(I) IS NOT CLOSE TO STRUE(I). 00163 C 00164 IF(.NOT. PASS) GO TO 5 00165 C PRINT FAIL MESSAGE AND HEADER. 00166 PASS = .FALSE. 00167 WRITE(NPRINT,1000) 00168 WRITE(NPRINT,1001) 00169 5 WRITE(NPRINT,1002) ICASE,N,INCX,INCY,MODE,I, 00170 * SCOMP(I),STRUE(I),SD,SSIZE(I) 00171 10 CONTINUE 00172 RETURN 00173 1000 FORMAT(1H+,39X,4HFAIL) 00174 1001 FORMAT(26H0CASE N INCX INCY MODE I, 00175 1 29X,7HCOMP(I),29X,7HTRUE(I),2X,10HDIFFERENCE, 00176 2 5X,7HSIZE(I)/1X) 00177 1002 FORMAT(1X,I4,I3,3I5,I3,2E36.8,2E12.4) 00178 END 00179 SUBROUTINE ITEST1(ICOMP,ITRUE) 00180 C1 ********************************* ITEST1 ************************* 00181 C 00182 C THIS SUBROUTINE COMPARES THE VARIABLES ICOMP AND ITRUE FOR 00183 C EQUALITY. 00184 C C. L. LAWSON, JPL, 1974 DEC 10 00185 C2 00186 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00187 LOGICAL PASS 00188 INTEGER ICOMP, ITRUE 00189 C 00190 IF(ICOMP .EQ. ITRUE) GO TO 10 00191 C 00192 C HERE ICOMP IS NOT EQUAL TO ITRUE. 00193 C 00194 IF(.NOT. PASS) GO TO 5 00195 C PRINT FAIL MESSAGE AND HEADER. 00196 PASS = .FALSE. 00197 WRITE(NPRINT,1000) 00198 WRITE(NPRINT,1001) 00199 5 ID=ICOMP-ITRUE 00200 WRITE(NPRINT,1002) ICASE,N,INCX,INCY,MODE,ICOMP,ITRUE,ID 00201 10 CONTINUE 00202 RETURN 00203 1000 FORMAT(1H+,39X,4HFAIL) 00204 1001 FORMAT(26H0CASE N INCX INCY MODE , 00205 1 29X,7HCOMP ,29X,7HTRUE ,2X,10HDIFFERENCE/1X) 00206 1002 FORMAT(1X,I4,I3,3I5,2I36,I12) 00207 END 00208 DOUBLE PRECISION FUNCTION DQDOTI(N,DB,QC,DX,INCX,DY,INCY) 00209 C D.P. DOT PRODUCT WITH EXTENDED PRECISION ACCUMULATION (AND RESULT) 00210 C QC AND DQDOTI ARE SET = DB + SUM FOR I = 0 TO N-1 OF 00211 C DX(LX+I*INCX) * DY(LY+I*INCY), WHERE QC IS AN EXTENDED 00212 C PRECISION RESULT WHICH CAN BE USED AS INPUT TO DQDOTA, 00213 C AND LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS 00214 C DEFINED IN A SIMILAR WAY USING INCY. THE MP PACKAGE BY 00215 C RICHARD P. BRENT IS USED FOR THE EXTENDED PRECISION ARITHMETIC. 00216 C 00217 C FRED T. KROGH, JPL, 1977, JUNE 1 00218 C2 00219 DOUBLE PRECISION DX(1), DY(1), DB 00220 INTEGER QC(10), QX(10), QY(10) 00221 C THE COMMON BLOCK FOR THE MP PACKAGE (MODIFIED TO GIVE IT A NAME) 00222 COMMON /MPCOM/ MPB, MPT, MPM, MPLUN, MPMXR, MPR(12) 00223 DATA I1 / 0 / 00224 C IF I1 IS 0 THE MP PACKAGE MUST BE INITIALIZED (MPBLAS SETS I1 = 1) 00225 IF (I1 .EQ. 0) CALL MPBLAS(I1) 00226 QC(1) = 0 00227 IF (DB .EQ. 0.D0) GO TO 60 00228 CALL MPCDM(DB, QX) 00229 CALL MPADD(QC, QX, QC) 00230 60 IF (N .EQ. 0) GO TO 80 00231 IX = 1 00232 IY = 1 00233 IF (INCX .LT. 0) IX = (-N + 1) * INCX + 1 00234 IF (INCY .LT. 0) IY = (-N + 1) * INCY + 1 00235 DO 70 I = 1,N 00236 CALL MPCDM(DX(IX), QX) 00237 CALL MPCDM(DY(IY), QY) 00238 CALL MPMUL(QX, QY, QX) 00239 CALL MPADD(QC, QX, QC) 00240 IX = IX + INCX 00241 IY = IY + INCY 00242 70 CONTINUE 00243 80 CALL MPCMD(QC, DQDOTI) 00244 RETURN 00245 END 00246 C PROGRAM TBLA 00247 C 00248 C THIS IS A TEST DRIVER FOR THE BLAS. 00249 C THE BLAS (BASIC LINEAR ALGEBRA SUBPROGRAMS) ARE A SET OF 00250 C THIRTY-EIGHT FORTRAN CALLABLE SUBPROGRAMS FOR BASIC OPERATIONS 00251 C OF NUMERICAL LINEAR ALGEBRA. THIS SOFTWARE PACKAGE IS THE 00252 C RESULT OF A VOLUNTARY AND COLLABORATIVE PROJECT OF THE 00253 C ACM-SIGNUM COMMITTEE ON BASIC LINEAR ALGEBRA SUBPROGRAMS. 00254 C THIS PROJECT WAS CARRIED OUT DURING THE PERIOD 1973-1977. 00255 C 00256 C THE BLAS ARE DESCRIBED IN THE PAPER, 00257 C BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, 00258 C BY C.L.LAWSON, R.J.HANSON, D.R.KINCAID, AND F.T.KROGH, 00259 C IN THE ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE,1979. 00260 C ALSO APPEARED AS U.TEXAS REPORT CNA-124, JULY, 1977, 00261 C AND SANDIA REPORT SAND77-0898J, FEBRUARY, 1978. 00262 C 00263 C 00264 C******************************************************************* 00265 C SUMMARY OF FUNCTIONS AND NAMES FOR BLAS 00266 C ------------------------------------------------------------------- 00267 C FUNCTION PREFIX AND SUFFIX ROOT 00268 C ------------------------------------------------------------------- 00269 C DOT PRODUCT /SDS- DS- DQ-I DQ-A C-U C-C D- S- -DOT 00270 C CONSTANT TIMES A VECTOR PLUS A VECTOR / C- D- S- -AXPY 00271 C SET-UP GIVENS ROTATION / D- S- -ROTG 00272 C APPLY ROTATION / D- S- -ROT 00273 C SET-UP MODIFIED GIVENS ROTATION / D- S- -ROTMG 00274 C APPLY MODIFIED ROTATION / D- S- -ROTM 00275 C COPY X TO Y / C- D- S- -COPY 00276 C SWAP X AND Y / C- D- S- -SWAP 00277 C 2-NORM (EUCLIDEAN LENGTH) / SC- D- S- -NRM2 00278 C SUM OF ABSOLUTE VALUES* / SC- D- S- -ASUM 00279 C CONSTANT TIMES A VECTOR / CS- C- D- S- -SCAL 00280 C INDEX OF ELEMENT HAVING MAX ABS VALUE*/ IC- ID- IS- -AMAX 00281 C ------------------------------------------------------------------ 00282 C *FOR COMPLEX VECTORS, THESE SUBPROGRAMS USE ABS(REAL)+ABS(IMAG). 00283 C 00284 C ARGUMENTS DESCRIBING VECTOR STORAGE 00285 C ----------------------------------- 00286 C 00287 C IN THE ARGUMENT LISTS, N DENOTES THE NUMBER OF COMPONENTS OF 00288 C A VECTOR, AND INCX DENOTES THE STORAGE SPACING BETWEEN COMPO- 00289 C NENTS OF THE X VECTOR. IF INCX .GE. 0 , THEN COMPONENT I OF 00290 C VECTOR X IS STORED IN SX(1+(I-1)*INCX) FOR I=1,...,N. 00291 C IF INCX .LT. 0 , COMPONENT I OF VECTOR X IS STORED IN 00292 C SX(1+(N-I)*IABS(INCX)). THE PARAMETER INCY GIVES THE STORAGE 00293 C SPACING FOR THE Y VECTOR. 00294 C ONLY POSITIVE VALUES OF INCX ARE ALLOWED FOR SUBPROGRAMS 00295 C THAT HAVE ONLY ONE VECTOR ARGUMENT. 00296 C 00297 C SPECIFICATION OF SUBPROGRAMS 00298 C ---------------------------- 00299 C DOT PRODUCT SUBPROGRAMS 00300 C ----------------------- 00301 C (SUM OF PRODUCTS OF COMPONENTS OF VECTORS X AND Y, 00302 C IF N .LE. 0 THE INNER PRODUCT WILL BE SET TO ZERO.) 00303 C SW = SDOT (N,SX,INCX,SY,INCY) 00304 C DW = DSDOT (N,SX,INCX,SY,INCY) 00305 C DOUBLE PRECISION ACCUMULATION USED IN DSDOT. 00306 C SW = SDSDOT (N,SB,SX,INCX,SY,INCY) 00307 C DOUBLE PRECISION ACCUMULATION AND DOUBLE PRECISION SUM OF 00308 C RESULTS PLUS SCALAR SB. SINGLE PRECISION RESULTS IN SW. 00309 C DW = DDOT (N,DX,INCX,DY,INCY) 00310 C DW = DQDOTI (N,DB,QC,DX,INCX,DY,INCY) 00311 C EXTENDED PRECISION ACCUMULATION AND EXTENDED PRECISION 00312 C SUM OF RESULTS PLUS DOUBLE PRECISION SCALAR DB. EXTENDED PRECISON 00313 C RESULTS IN QC AND DOUBLE PRECISION RESULTS IN DW. 00314 C DW = DQDOTA (N,DB,QC,DX,INCX,DY,INCY) 00315 C EXTENDED PRECISION ACCUMULATION AND EXTENDED PRECISION 00316 C SUM OF RESULTS PLUS EXTENDED PRECISION SCALAR QC AND DOUBLE PRECISION 00317 C SCALAR DB. EXTENDED PRECISION RESULTS IN QC AND DOUBLE 00318 C PRECISON RESULTS IN DW. 00319 C CW = CDOTC (N,CX,INCX,CY,INCY) 00320 C COMPLEX CONJUGATE OF X VECTOR USED. 00321 C CW = CDOTU (N,CX,INCX,CY,INCY) 00322 C UNCONJUGATED VECTORS USED. 00323 C 00324 C ELEMENTARY VECTOR OPERATION (Y = A*X + Y) 00325 C ----------------------------------------- 00326 C CALL SAXPY (N,SA,SX,INCX,SY,INCY) 00327 C CALL DAXPY (N,DA,DX,INCX,DY,INCY) 00328 C CALL CAXPY (N,CA,CX,INCX,CY,INCY) 00329 C IF A=0 OR IF N .LE. 0 THESE SUBROUTINES RETURN IMMEDIATELY. 00330 C 00331 C CONSTRUCT GIVENS PLANE ROTATION 00332 C ------------------------------- 00333 C CALL SROTG (SA,SB,SC,SS) 00334 C CALL DROTG (DA,DB,DC,DS) 00335 C SEE TOMS PAPER FOR DETAILS. 00336 C 00337 C APPLY A PLANE ROTATION 00338 C ---------------------- 00339 C CALL SROT (N,SX,INCX,SY,INCY,SC,SS) 00340 C CALL DROT (N,DX,INCX,DY,INCY,DC,DS) 00341 C SEE TOMS PAPER FOR DETAILS. 00342 C 00343 C CONSTRUCT A MODIFIED GIVENS TRANSFORMATION 00344 C ------------------------------------------ 00345 C CALL SROTMG (SD1,SD2,SB1,SB2,SPARAM) 00346 C CALL DROTMG (DD1,DD2,DB1,DB2,DPARAM) 00347 C SEE TOMS PAPER FOR DETAILS. 00348 C 00349 C APPLY A MODIFIED GIVENS TRANSFORMATION 00350 C -------------------------------------- 00351 C CALL SROTM (N,SX,INCX,SY,INCY,SPARAM) 00352 C CALL DROTM (N,DX,INCX,DY,INCY,DPARAM) 00353 C SEE TOMS PAPER FOR DETAILS. 00354 C 00355 C COPY A VECTOR X TO Y 00356 C -------------------- 00357 C CALL SCOPY (N,SX,INCX,SY,INCY) 00358 C CALL DCOPY (N,DX,INCX,DY,INCY) 00359 C CALL CCOPY (N,CX,INCX,CY,INCY) 00360 C IF N .LE. 0 THESE SUBROUTINES RETURN IMMEDIATELY 00361 C 00362 C INTERCHANGE VECTORS X AND Y 00363 C --------------------------- 00364 C CALL SSWAP (N,SX,INCX,SY,INCY) 00365 C CALL DSWAP (N,DX,INCX,DY,INCY) 00366 C CALL CSWAP (N,CX,INCX,CY,INCY) 00367 C IF N .LE. 0 THESE SUBROUTINES RETURN IMMEDIATELY 00368 C 00369 C EUCLIDEAN LENGTH OR L-2 NORM OF A VECTOR 00370 C ---------------------------------------- 00371 C (SQUARE ROOT OF SUM OF ABSOLUTE VALUES SQUARED.) 00372 C SW = SNRM2 (N,SX,INCX) 00373 C DW = DNRM2 (N,DX,INCX) 00374 C SW = SCNRM2 (N,CX,INCX) 00375 C IF N .LE. THESE SUBROUTINES RETURN IMMEDIATELY 00376 C 00377 C SUM OF MAGNITUDES OF VECTOR COMPONENTS 00378 C -------------------------------------- 00379 C (SUM OF ABSOLUTE VALUES OR ABS(REAL)+ABS(IMAG)) 00380 C SW = SASUM (N,SX,INCX) 00381 C DW = DASUM (N,DX,INCX) 00382 C SW = SCASUM (N,CX,INCX) 00383 C IF N .LE. 0 THESE FUNCTIONS ARE SET TO 0 AND RETURN IMMEDIATELY. 00384 C 00385 C VECTOR SCALING (X = A*X) 00386 C ------------------------- 00387 C CALL SSCAL (N,SA,SX,INCX) 00388 C CALL DSCAL (N,DA,DX,INCX) 00389 C CALL CSCAL (N,CA,CX,INCX) 00390 C CALL CSSCAL (N,SA,CX,INCX) 00391 C IF N .LE. 0 THESE SUBPROGRAMS RETURN IMMEDIATELY. 00392 C 00393 C FIND LARGEST COMPONENT OF A VECTOR 00394 C ---------------------------------- 00395 C (SMALLEST INDEX OF COMPONENT WITH LARGEST ABSOLUTE VALUE OR 00396 C ABS(REAL)+ABS(IMAG).) 00397 C IMAX = ISAMAX (N,SX,INCX) 00398 C IMAX = IDAMAX (N,DX,INCX) 00399 C IMAX = ICAMAX (N,CX,INCX) 00400 C IF N .LE. 0 THESE FUNCTIONS SET TO 0 AND RETURN IMMEDIATELY. 00401 C 00402 C TYPE DECLARATIONS FOR FUNCTION NAMES ARE AS FOLLOWS.. 00403 C 00404 C INTEGER ISAMAX,IDAMAX,ICAMAX 00405 C REAL SDOT,SDSDOT,SNRM2,SCNRM2,SASUM,SCASUM 00406 C DOUBLE PRECISION DSDOT,DDOT,DQDOTI,DQDOTA,DASUM 00407 C COMPLEX CDOTC,CDOTU 00408 C 00409 C TYPE AND DIMENSION INFORMATION FOR VARIABLES OCCURRING IN 00410 C SUBPROGRAM SPECIFICATIONS ARE AS FOLLOWS.. 00411 C 00412 C INTEGER N,INXC,INCY,IMAX 00413 C REAL SC(MX),SY(MY),SA,SB,SC,SS 00414 C REAL SD1,SD2,SB1,SB2,SPARAM(5),SW,QC(10) 00415 C DOUBLE PRECISION DX(MX),DY(MY),DA,DB,DC,DS 00416 C DOUBLE PRECISION DD1,DD2,DB1,DB2,DPARAM(5),DW 00417 C COMPLEX CX(MX),CY(MY),CA,CW 00418 C 00419 C WHERE MX = MAX(1,N*ABS(INCX)) 00420 C MY = MAX(1,N*ABS(INCY)) 00421 C 00422 C 00423 C************* DEMONSTRATION OF USAGE OF BLAS ********************** 00424 C 00425 C DIMENSION A(20,20),B(15,10),C(20,15),X(10) 00426 C INTEGER IP(20) 00427 C 00428 C MDA = 20 00429 C MDB = 15 00430 C MDC = 20 00431 C 00432 C M = 10 00433 C K = 15 00434 C N = 10 00435 C 00436 C------------------------------------------------------------------- 00437 C PRODUCT OF RECTANGULAR MATRICES C(MXN) = A(MXK)*B(KXN) 00438 C 00439 C DO 10 J=1,M 00440 C DO 10 I=1,N 00441 C 10 C(I,J) = SDOT(K,A(I,1),MDA,B(1,J),1) 00442 C 00443 C------------------------------------------------------------------- 00444 C SOLVE NXN UPPER TRANGULAR NONSINGULAR LINEAR SYSTEM AX = B 00445 C 00446 C DO 20 II=1,N 00447 C I = N+1-II 00448 C CALL SSCAL(M,1./A(I,I),B(I,1),MDB) 00449 C DO 20 J=1,M 00450 C 20 CALL SAXPY(I-1,-B(I,J),A(1,I),1,B(1,J),1) 00451 C 00452 C------------------------------------------------------------------- 00453 C SCALE COLUMNS OF RECTANGULAR MATRIX C(MXN) 00454 C 00455 C DO 30 J=1,N 00456 C T = 1.E0/SNRM2(M,C(1,J),1) 00457 C 30 CALL SSCAL(M,T,C(1,J),1) 00458 C 00459 C------------------------------------------------------------------- 00460 C ROW EQUILIBRATE SQUARE MATRIX A(NXN) 00461 C 00462 C DO 40 I=1,N 00463 C JMAX = ISAMAX(N,A(I,1),MDA) 00464 C T = A(I,JMAX) 00465 C IF(T .EQ. 0.E0) GO TO 40 00466 C CALL SSCAL(N,1.E0/T,A(I,1),MDA) 00467 C 40 CONTINUE 00468 C 00469 C----------------------------------------------------------------- 00470 C TO CHOOSE ROW PIVOT IN GAUSSIAN ELIMINATION USE 00471 C 00472 C IMAX = ISAMAX(N-J+1,A(J,J),1) + J-1 00473 C 00474 C------------------------------------------------------------------- 00475 C SET NXN MATRIX A TO IDENTITY MATRIX AND SET B = A 00476 C 00477 C DO 50 J=1,N 00478 C 50 CALL SCOPY(N,0.E0,0,A(1,J),1) 00479 C CALL SCOPY(N,1.E0,0,A,MDA+1) 00480 C DO 60 J=1,N 00481 C 60 CALL SCOPY(N,A(1,J),1,B(1,J),1) 00482 C 00483 C------------------------------------------------------------------- 00484 C INTERCHANGE OR SWAP COLUMNS OF MXN MATRIX C 00485 C 00486 C DO 70 J=1,N 00487 C L = IP(J) 00488 C IF(J .NE. L) CALL SSWAP(M,C(1,J),1,C(1,L),1) 00489 C 70 CONTINUE 00490 C 00491 C------------------------------------------------------------------- 00492 C TRANSPOSE AN NXN MATRIX A IN-PLACE 00493 C 00494 C DO 80 J=1,N 00495 C 80 CALL SSWAP(N-J,A(J,J+1),MDA,A(J+1,J),1) 00496 C 00497 C END 00498 C PROGRAM MAIN(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) 00499 C1 ********************************* TBLA *************************** 00500 C TEST DRIVER FOR BASIC LINEAR ALGEBRA SUBPROGRAMS. 00501 C C. L. LAWSON,JPL, 1974 DEC 10, 1975 MAY 28 00502 C TBLA READS INPUT FROM UNIT 5 AND WRITES OUTPUT ON UNIT 00503 C NPRINT WHICH IS NOMINALLY SET TO 6. THE FORM OF ACCEPTABLE 00504 C INPUT IS DESCRIBED IN FORMAT STATEMENT 1002. 00505 C FOR EACH SUBPROGRAM SELECTED FOR TESTING, TBLA CALLS ONE OF 00506 C THE SUBROUTINES CHECK0, CHECK1, CHECK2. CHECK0 IS USED TO TEST 00507 C SUBPROGRAMS HAVING NO VECTOR ARGUMENTS, CHECK1 FOR THOSE HAVING 00508 C ONE VECTOR ARGUMENT, AND CHECK2 FOR THOSE HAVING TWO. 00509 C THE TUNING PARAMETERS SFAC, SDFAC, DFAC, AND DQFAC ARE SET IN 00510 C A DATA STATEMENT AND PASSED TO CHECK0, CHECK1, AND CHECK2 TO SET 00511 C TOLERANCES ON TESTING THE SUBPROGRAMS. THE PREFIXES, S,SD,D, 00512 C AND DQ REFER TO THE TYPE OF SUBPROGRAM FOR WHICH EACH TOLERANCE 00513 C IS USED, NAMELY SINGLE PRECISION, MIXED SINGLE AND DOUBLE 00514 C PRECISION, DOUBLE PRECISION, AND MIXED DOUBLE AND EXTENDED 00515 C PRECISION. 00516 C THE TUNING PARAMETERS ULTIMATELY ARE USED IN STEST AND DTEST. 00517 C SEE THESE SUBROUTINE LISTINGS FOR THE PRECISE ROLE OF THOSE 00518 C PARAMETERS. THESE PARAMETERS COMPENSATE FOR THE VAGARIES OF 00519 C ARITHMETIC TRUNCATION ON DIFFERENT MACHINES. SETTING A TUNING 00520 C PARAMETER SMALLER PROVIDES MORE TOLERANCE FOR BAD TRUNCATION, 00521 C I.E. MAKES IT EASIER TO PASS THE TESTS. 00522 C THE PARAMETERS IN COMMON/COMBLA/ ARE USED AS FOLLOWS.. 00523 C 00524 C NPRINT FORTRAN UNIT FOR PRINTED OUTPUT. SET IN TBLA. 00525 C USED IN TBLA, HEADER, STEST, DTEST, AND ITEST1. 00526 C ICASE NUMBER IDENTIFYING SUBPROGRAM BEING TESTED. SEE COMMENTS 00527 C ALONG RIGHT MARGIN IN CHECK0, CHECK1, AND CHECK2 00528 C FOR ASSOCIATION OF NUMBERS FROM 1 TO 38 WITH NAMES OF 00529 C SUBPROGRAMS. ICASE IS SET IN TBLA AND USED IN VARIOUS 00530 C OF THE SUBROUTINES. 00531 C N SET IN CHECK0, CHECK1, OR CHECK2. GENERALLY DENOTES 00532 C THE DIMENSION OF A VECTOR BEING SENT TO A BLAS 00533 C SUBPROGRAM, BUT IN TESTS NOT INVOLVING VECTOR 00534 C ARGUMENTS N IS USED JUST TO DISTINGUISH DIFFERENT SETS 00535 C OF TEST DATA. WILL BE PRINTED WHEN ERRORS ARE NOTED. 00536 C INCX SET IN TBLA, CHECK1, AND CHECK2. SENT TO BLAS 00537 C SUBPROGRAMS AS TEST DATA. PRINTED WHEN ERRORS ARE 00538 C NOTED. 00539 C INCY SET IN TBLA, AND CHECK2. SENT TO BLAS SUBPROGRAMS AS 00540 C TEST DATA. PRINTED WHEN ERRORS ARE NOTED. 00541 C MODE SET IN TBLA AND CHECK2. DISTINGUISHES TEST CASES. 00542 C PRINTED WHEN ERRORS ARE NOTED. 00543 C PASS SET IN TBLA, STEST, DTEST, AND ITEST1. SET TO TRUE 00544 C OR FALSE TO DENOTE SUCCESS OR FAILURE OF TESTING FOR 00545 C A BLAS SUBPROGRAM. ALWAYS PRINTED FOR EACH SUBPROGRAM 00546 C TESTED. 00547 C2 00548 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00549 LOGICAL PASS,ALLZRO 00550 INTEGER ITEST(38) 00551 DOUBLE PRECISION DFAC,DQFAC 00552 DATA SFAC,SDFAC,DFAC,DQFAC / .3125E-1, .50, .625D-1, .125D0/ 00553 NPRINT = 6 00554 5 WRITE(NPRINT,1002) 00555 READ(5,1000) ITEST 00556 WRITE(NPRINT,1005) ITEST 00557 ALLZRO=.TRUE. 00558 DO 60 IC=1,38 00559 ICASE=IC 00560 IF(ITEST(ICASE) .EQ. 0) GO TO 60 00561 ALLZRO=.FALSE. 00562 CALL HEADER 00563 C 00564 C INITIALIZE PASS, INCX, INCY, AND MODE FOR A NEW CASE. 00565 C THE VALUE 9999 FOR INCX, INCY OR MODE WILL APPEAR IN THE 00566 C DETAILED OUTPUT, IF ANY, FOR CASES THAT DO NOT INVOLVE 00567 C THESE PARAMETERS. 00568 C 00569 PASS=.TRUE. 00570 INCX=9999 00571 INCY=9999 00572 MODE=9999 00573 GO TO (12,12,12,12,12,12,12,12,12,12, 00574 A 12,10,10,12,12,10,10,12,12,12, 00575 B 12,12,12,12,12,11,11,11,11,11, 00576 C 11,11,11,11,11,11,11,11), ICASE 00577 C ICASE = 12-13 OR 16-17 00578 10 CALL CHECK0(SFAC,DFAC) 00579 GO TO 50 00580 C ICASE = 26-38 00581 11 CALL CHECK1(SFAC,DFAC) 00582 GO TO 50 00583 C ICASE = 1-11, 14-15, OR 18-25 00584 12 CALL CHECK2(SFAC,SDFAC,DFAC,DQFAC) 00585 50 CONTINUE 00586 C PRINT 00587 IF(PASS) WRITE(NPRINT,1001) 00588 60 CONTINUE 00589 IF(.NOT. ALLZRO) GO TO 5 00590 STOP 00591 1000 FORMAT(80I1) 00592 1001 FORMAT(1H+,39X,4HPASS) 00593 1002 FORMAT(1X///33H PROGRAM TBLA IS READY FOR INPUT./ 00594 142H INPUT ONE CARD IMAGE HAVING ONES OR ZEROS/ 00595 242H IN COLS 1 - 38. A ONE IN COL K MEANS TO/ 00596 341H TEST SUBPROGRAM NO. K. ALL ZEROS MEANS/ 00597 438H TO TERMINATE EXECUTION. INPUT NOW.) 00598 1005 FORMAT(1H0,38I2) 00599 END 00600 SUBROUTINE HEADER 00601 C1 ********************************* HEADER ************************* 00602 C PRINT HEADER FOR CASE 00603 C C. L. LAWSON, JPL, 1974 DEC 12 00604 C2 00605 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00606 LOGICAL PASS 00607 DIMENSION L(3,38) 00608 C 00609 DATA L(1, 1),L(2, 1),L(3, 1)/2H ,2HSD,2HOT/ 00610 DATA L(1, 2),L(2, 2),L(3, 2)/2H D,2HSD,2HOT/ 00611 DATA L(1, 3),L(2, 3),L(3, 3)/2HSD,2HSD,2HOT/ 00612 DATA L(1, 4),L(2, 4),L(3, 4)/2H ,2HDD,2HOT/ 00613 DATA L(1, 5),L(2, 5),L(3, 5)/2HDQ,2HDO,2HTI/ 00614 DATA L(1, 6),L(2, 6),L(3, 6)/2HDQ,2HDO,2HTA/ 00615 DATA L(1,7),L(2,7),L(3,7)/2H C,2HDO,2HTC/ 00616 DATA L(1, 8),L(2, 8),L(3, 8)/2H C,2HDO,2HTU/ 00617 DATA L(1, 9),L(2, 9),L(3, 9)/2H S,2HAX,2HPY/ 00618 DATA L(1,10),L(2,10),L(3,10)/2H D,2HAX,2HPY/ 00619 DATA L(1,11),L(2,11),L(3,11)/2H C,2HAX,2HPY/ 00620 DATA L(1,12),L(2,12),L(3,12)/2H S,2HRO,2HTG/ 00621 DATA L(1,13),L(2,13),L(3,13)/2H D,2HRO,2HTG/ 00622 DATA L(1,14),L(2,14),L(3,14)/2H ,2HSR,2HOT/ 00623 DATA L(1,15),L(2,15),L(3,15)/2H ,2HDR,2HOT/ 00624 DATA L(1,16),L(2,16),L(3,16)/2HSR,2HOT,2HMG/ 00625 DATA L(1,17),L(2,17),L(3,17)/2HDR,2HOT,2HMG/ 00626 DATA L(1,18),L(2,18),L(3,18)/2H S,2HRO,2HTM/ 00627 DATA L(1,19),L(2,19),L(3,19)/2H D,2HRO,2HTM/ 00628 DATA L(1,20),L(2,20),L(3,20)/2H S,2HCO,2HPY/ 00629 DATA L(1,21),L(2,21),L(3,21)/2H D,2HCO,2HPY/ 00630 DATA L(1,22),L(2,22),L(3,22)/2H C,2HCO,2HPY/ 00631 DATA L(1,23),L(2,23),L(3,23)/2H S,2HSW,2HAP/ 00632 DATA L(1,24),L(2,24),L(3,24)/2H D,2HSW,2HAP/ 00633 DATA L(1,25),L(2,25),L(3,25)/2H C,2HSW,2HAP/ 00634 DATA L(1,26),L(2,26),L(3,26)/2H S,2HNR,2HM2/ 00635 DATA L(1,27),L(2,27),L(3,27)/2H D,2HNR,2HM2/ 00636 DATA L(1,28),L(2,28),L(3,28)/2HSC,2HNR,2HM2/ 00637 DATA L(1,29),L(2,29),L(3,29)/2H S,2HAS,2HUM/ 00638 DATA L(1,30),L(2,30),L(3,30)/2H D,2HAS,2HUM/ 00639 DATA L(1,31),L(2,31),L(3,31)/2HSC,2HAS,2HUM/ 00640 DATA L(1,32),L(2,32),L(3,32)/2H S,2HSC,2HAL/ 00641 DATA L(1,33),L(2,33),L(3,33)/2H D,2HSC,2HAL/ 00642 DATA L(1,34),L(2,34),L(3,34)/2H C,2HSC,2HAL/ 00643 DATA L(1,35),L(2,35),L(3,35)/2HCS,2HSC,2HAL/ 00644 DATA L(1,36),L(2,36),L(3,36)/2HIS,2HAM,2HAX/ 00645 DATA L(1,37),L(2,37),L(3,37)/2HID,2HAM,2HAX/ 00646 DATA L(1,38),L(2,38),L(3,38)/2HIC,2HAM,2HAX/ 00647 C 00648 WRITE(NPRINT,1000) ICASE,(L(I,ICASE),I=1,3) 00649 RETURN 00650 C 00651 1000 FORMAT(23H0TEST OF SUBPROGRAM NO.,I3,2X,3A2) 00652 END 00653 SUBROUTINE CHECK0(SFAC,DFAC) 00654 C1 ********************************* CHECK0 ************************* 00655 C THIS SUBROUTINE TESTS SUBPROGRAMS 12-13 AND 16-17. 00656 C THESE SUBPROGRAMS HAVE NO ARRAY ARGUMENTS. 00657 C 00658 C C. L. LAWSON, JPL, 1975 MAR 07, MAY 28 00659 C R. J. HANSON, J. A. WISNIEWSKI, SANDIA LABS, APRIL 25,1977. 00660 C2 00661 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00662 LOGICAL PASS 00663 REAL STRUE(9),STEMP(9) 00664 DOUBLE PRECISION DC,DS,DA1(8),DB1(8),DC1(8),DS1(8) 00665 DOUBLE PRECISION DA,DATRUE(8),DBTRUE(8),DZERO,DFAC,DB 00666 DOUBLE PRECISION DAB(4,9),DTEMP(9),DTRUE(9,9),D12 00667 DATA ZERO, DZERO / 0., 0.D0 / 00668 DATA DA1/ .3D0, .4D0, -.3D0, -.4D0, -.3D0, 0.D0, 0.D0, 1.D0/ 00669 DATA DB1/ .4D0, .3D0, .4D0, .3D0, -.4D0, 0.D0, 1.D0, 0.D0/ 00670 DATA DC1/ .6D0, .8D0, -.6D0, .8D0, .6D0, 1.D0, 0.D0, 1.D0/ 00671 DATA DS1/ .8D0, .6D0, .8D0, -.6D0, .8D0, 0.D0, 1.D0, 0.D0/ 00672 DATA DATRUE/ .5D0, .5D0, .5D0, -.5D0, -.5D0, 0.D0, 1.D0, 1.D0/ 00673 DATA DBTRUE/ 0.D0, .6D0, 0.D0, -.6D0, 0.D0, 0.D0, 1.D0, 0.D0/ 00674 C INPUT FOR MODIFIED GIVENS 00675 DATA DAB/ .1D0,.3D0,1.2D0,.2D0, 00676 A .7D0, .2D0, .6D0, 4.2D0, 00677 B 0.D0,0.D0,0.D0,0.D0, 00678 C 4.D0, -1.D0, 2.D0, 4.D0, 00679 D 6.D-10, 2.D-2, 1.D5, 10.D0, 00680 E 4.D10, 2.D-2, 1.D-5, 10.D0, 00681 F 2.D-10, 4.D-2, 1.D5, 10.D0, 00682 G 2.D10, 4.D-2, 1.D-5, 10.D0, 00683 H 4.D0, -2.D0, 8.D0, 4.D0 / 00684 C TRUE RESULTS FOR MODIFIED GIVENS 00685 DATA DTRUE/0.D0,0.D0, 1.3D0, .2D0, 0.D0,0.D0,0.D0, .5D0, 0.D0, 00686 A 0.D0,0.D0, 4.5D0, 4.2D0, 1.D0, .5D0, 0.D0,0.D0,0.D0, 00687 B 0.D0,0.D0,0.D0,0.D0, -2.D0, 0.D0,0.D0,0.D0,0.D0, 00688 C 0.D0,0.D0,0.D0, 4.D0, -1.D0, 0.D0,0.D0,0.D0,0.D0, 00689 D 0.D0, 15.D-3, 0.D0, 10.D0, -1.D0, 0.D0, -1.D-4, 00690 E 0.D0, 1.D0, 00691 F 0.D0,0.D0, 6144.D-5, 10.D0, -1.D0, 4096.D0, -1.D6, 00692 G 0.D0, 1.D0, 00693 H 0.D0,0.D0,15.D0,10.D0,-1.D0, 5.D-5, 0.D0,1.D0,0.D0, 00694 I 0.D0,0.D0, 15.D0, 10.D0, -1. D0, 5.D5, -4096.D0, 00695 J 1.D0, 4096.D-6, 00696 K 0.D0,0.D0, 7.D0, 4.D0, 0.D0,0.D0, -.5D0, -.25D0, 0.D0/ 00697 C 4096 = 2 ** 12 00698 DATA D12 /4096.D0/ 00699 C 00700 C COMPUTE TRUE VALUES WHICH CANNOT BE PRESTORED 00701 C IN DECIMAL NOTATION. 00702 DTRUE(1,1) = 12.D0 / 130.D0 00703 DTRUE(2,1) = 36.D0 / 130.D0 00704 DTRUE(7,1) = -1.D0 / 6.D0 00705 DTRUE(1,2) = 14.D0 / 75.D0 00706 DTRUE(2,2) = 49.D0 / 75.D0 00707 DTRUE(9,2) = 1.D0 / 7.D0 00708 DTRUE(1,5) = 45.D-11 * (D12 * D12) 00709 DTRUE(3,5) = 4.D5 / (3.D0 * D12) 00710 DTRUE(6,5) = 1.D0 / D12 00711 DTRUE(8,5) = 1.D4 / (3.D0 * D12) 00712 DTRUE(1,6) = 4.D10 / (1.5D0 * D12 * D12) 00713 DTRUE(2,6) = 2.D-2 / 1.5D0 00714 DTRUE(8,6) = 5.D-7 * D12 00715 DTRUE(1,7) = 4.D0 / 150.D0 00716 DTRUE(2,7) = (2.D-10 / 1.5D0) * (D12 * D12) 00717 DTRUE(7,7) = -DTRUE(6,5) 00718 DTRUE(9,7) = 1.D4 / D12 00719 DTRUE(1,8) = DTRUE(1,7) 00720 DTRUE(2,8) = 2.D10 / (1.5D0 * D12 * D12) 00721 DTRUE(1,9) = 32.D0 / 7.D0 00722 DTRUE(2,9) = -16.D0 / 7.D0 00723 DBTRUE(1) = 1.D0/.6D0 00724 DBTRUE(3) = -1.D0/.6D0 00725 DBTRUE(5) = 1.D0/.6D0 00726 C 00727 JUMP= ICASE-11 00728 DO 500 K=1,9 00729 C SET N=K FOR IDENTIFICATION IN OUTPUT IF ANY. 00730 N=K 00731 C BRANCH TO SELECT SUBPROGRAM TO BE TESTED. 00732 C 00733 GO TO (120,130,999,999,160,170), JUMP 00734 C 12. SROTG 00735 120 IF(K.GT.8) GO TO 600 00736 SA=SNGL(DA1(K)) 00737 SB = SNGL(DB1(K)) 00738 CALL SROTG(SA,SB,SC,SS) 00739 CALL STEST1 (SA,SNGL(DATRUE(K)),SNGL(DATRUE(K)),SFAC) 00740 CALL STEST1 (SB,SNGL(DBTRUE(K)),SNGL(DBTRUE(K)),SFAC) 00741 CALL STEST1 (SC,SNGL(DC1(K)),SNGL(DC1(K)),SFAC) 00742 CALL STEST1 (SS,SNGL(DS1(K)),SNGL(DS1(K)),SFAC) 00743 GO TO 500 00744 C 13. DROTG 00745 130 IF(K.GT.8) GO TO 600 00746 DA=DA1(K) 00747 DB = DB1(K) 00748 CALL DROTG(DA,DB,DC,DS) 00749 CALL DTEST1 (DA,DATRUE(K),DATRUE(K),DFAC) 00750 CALL DTEST1 (DB,DBTRUE(K),DBTRUE(K),DFAC) 00751 CALL DTEST1 (DC,DC1(K),DC1(K),DFAC) 00752 CALL DTEST1 (DS,DS1(K),DS1(K),DFAC) 00753 GO TO 500 00754 C 16. SROTMG 00755 160 CONTINUE 00756 DO 162 I=1,4 00757 STEMP(I) = SNGL(DAB(I,K)) 00758 STEMP(I+4) = ZERO 00759 162 CONTINUE 00760 STEMP(9) = ZERO 00761 CALL SROTMG(STEMP(1),STEMP(2),STEMP(3),STEMP(4),STEMP(5)) 00762 C 00763 DO 166 I=1,9 00764 166 STRUE(I)=SNGL(DTRUE(I,K)) 00765 CALL STEST(9,STEMP,STRUE,STRUE,SFAC) 00766 GO TO 500 00767 C 17. DROTMG 00768 170 CONTINUE 00769 DO 172 I=1,4 00770 DTEMP(I)= DAB(I,K) 00771 DTEMP(I+4) = DZERO 00772 172 CONTINUE 00773 DTEMP(9) = DZERO 00774 CALL DROTMG(DTEMP(1),DTEMP(2),DTEMP(3),DTEMP(4),DTEMP(5)) 00775 DO I=1,9 00776 IF (ABS(DTEMP(I)-DTRUE(I,K)) .GT. 0.0001 ) THEN 00777 WRITE(*,*) "DTEMP-DTRUE=",I,DTEMP(I)-DTRUE(I,K) 00778 END IF 00779 END DO 00780 C CALL DTEST(9,DTEMP,DTRUE(1,K),DTRUE(1,K),DFAC) 00781 500 CONTINUE 00782 600 RETURN 00783 C THE FOLLOWING STOP SHOULD NEVER BE REACHED. 00784 999 STOP 00785 END 00786 SUBROUTINE CHECK1(SFAC,DFAC) 00787 C1 ********************************* CHECK1 ************************* 00788 C THIS SUBPROGRAM TESTS THE INCREMENTING AND ACCURACY OF THE LINEAR 00789 C ALGEBRA SUBPROGRAMS 26 - 38 (SNRM2 TO ICAMAX). STORED RESULTS ARE 00790 C COMPARED WITH THE RESULT RETURNED BY THE SUBPROGRAM. 00791 C 00792 C THESE SUBPROGRAMS REQUIRE A SINGLE VECTOR ARGUMENT. 00793 C 00794 C ICASE DESIGNATES WHICH SUBPROGRAM TO TEST. 00795 C 26 .LE. ICASE .LE. 38 00796 C C. L. LAWSON, JPL, 1974 DEC 10, MAY 28 00797 C2 00798 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00799 LOGICAL PASS 00800 INTEGER ITRUE2(5),ITRUE3(5) 00801 DOUBLE PRECISION DA,DX(8) 00802 DOUBLE PRECISION DV(8,5,2) 00803 DOUBLE PRECISION DFAC 00804 DOUBLE PRECISION DNRM2,DASUM 00805 DOUBLE PRECISION DTRUE1(5),DTRUE3(5),DTRUE5(8,5,2) 00806 REAL STRUE2(5),STRUE4(5),STRUE(8),SX(8) 00807 COMPLEX CA,CV(8,5,2),CTRUE5(8,5,2),CTRUE6(8,5,2),CX(8) 00808 C 00809 DATA SA, DA, CA / .3, .3D0, (.4,-.7) / 00810 DATA DV/.1D0,2.D0,2.D0,2.D0,2.D0,2.D0,2.D0,2.D0, 00811 1 .3D0,3.D0,3.D0,3.D0,3.D0,3.D0,3.D0,3.D0, 00812 2 .3D0,-.4D0,4.D0,4.D0,4.D0,4.D0,4.D0,4.D0, 00813 3 .2D0,-.6D0,.3D0,5.D0,5.D0,5.D0,5.D0,5.D0, 00814 4 .1D0,-.3D0,.5D0,-.1D0,6.D0,6.D0,6.D0,6.D0, 00815 5 .1D0,8.D0,8.D0,8.D0,8.D0,8.D0,8.D0,8.D0, 00816 6 .3D0,9.D0,9.D0,9.D0,9.D0,9.D0,9.D0,9.D0, 00817 7 .3D0,2.D0,-.4D0,2.D0,2.D0,2.D0,2.D0,2.D0, 00818 8 .2D0,3.D0,-.6D0,5.D0,.3D0,2.D0,2.D0,2.D0, 00819 9 .1D0,4.D0,-.3D0,6.D0,-.5D0,7.D0,-.1D0, 3.D0/ 00820 C COMPLEX TEST VECTORS 00821 DATA CV/ 00822 1(.1,.1),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.), 00823 2(.3,-.4),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.), 00824 3(.1,-.3),(.5,-.1),(5.,6.),(5.,6.),(5.,6.),(5.,6.),(5.,6.),(5.,6.), 00825 4(.1,.1),(-.6,.1),(.1,-.3),(7.,8.),(7.,8.),(7.,8.),(7.,8.),(7.,8.), 00826 5(.3,.1),(.1,.4),(.4,.1),(.1,.2),(2.,3.),(2.,3.),(2.,3.),(2.,3.), 00827 6(.1,.1),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.), 00828 7(.3,-.4),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.), 00829 8(.1,-.3),(8.,9.),(.5,-.1),(2.,5.),(2.,5.),(2.,5.),(2.,5.),(2.,5.), 00830 9(.1,.1),(3.,6.),(-.6,.1),(4.,7.),(.1,-.3),(7.,2.),(7.,2.),(7.,2.), 00831 T(.3,.1),(5.,8.),(.1,.4),(6.,9.),(.4,.1),(8.,3.),(.1,.2),(9.,4.) / 00832 C 00833 DATA STRUE2/.0,.5,.6,.7,.7/ 00834 DATA STRUE4/.0,.7,1.,1.3,1.7/ 00835 DATA DTRUE1/.0D0,.3D0,.5D0,.7D0,.6D0/ 00836 DATA DTRUE3/.0D0,.3D0,.7D0,1.1D0,1.D0/ 00837 DATA DTRUE5/.10D0,2.D0,2.D0,2.D0,2.D0,2.D0,2.D0,2.D0, 00838 1 .09D0,3.D0,3.D0,3.D0,3.D0,3.D0,3.D0,3.D0, 00839 2 .09D0,-.12D0,4.D0,4.D0,4.D0,4.D0,4.D0,4.D0, 00840 3 .06D0,-.18D0,.09D0,5.D0,5.D0,5.D0,5.D0,5.D0, 00841 4 .03D0,-.09D0,.15D0,-.03D0,6.D0,6.D0,6.D0,6.D0, 00842 5 .10D0,8.D0,8.D0,8.D0,8.D0,8.D0,8.D0,8.D0, 00843 6 .09D0,9.D0,9.D0,9.D0,9.D0,9.D0,9.D0,9.D0, 00844 7 .09D0,2.D0,-.12D0,2.D0,2.D0,2.D0,2.D0,2.D0, 00845 8 .06D0,3.D0,-.18D0,5.D0,.09D0,2.D0,2.D0,2.D0, 00846 9 .03D0,4.D0, -.09D0,6.D0, -.15D0,7.D0, -.03D0, 3.D0/ 00847 C 00848 DATA CTRUE5/ 00849 A(.1,.1),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.), 00850 B(-.16,-.37),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.), 00851 C (3.,4.), 00852 D(-.17,-.19),(.13,-.39),(5.,6.),(5.,6.),(5.,6.),(5.,6.),(5.,6.), 00853 E (5.,6.), 00854 F(.11,-.03),(-.17,.46),(-.17,-.19),(7.,8.),(7.,8.),(7.,8.),(7.,8.), 00855 G (7.,8.), 00856 H(.19,-.17),(.32,.09),(.23,-.24),(.18,.01),(2.,3.),(2.,3.),(2.,3.), 00857 I (2.,3.), 00858 J(.1,.1),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.), 00859 K(-.16,-.37),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.), 00860 L (6.,7.), 00861 M(-.17,-.19),(8.,9.),(.13,-.39),(2.,5.),(2.,5.),(2.,5.),(2.,5.), 00862 N (2.,5.), 00863 O(.11,-.03),(3.,6.),(-.17,.46),(4.,7.),(-.17,-.19),(7.,2.),(7.,2.), 00864 P (7.,2.), 00865 Q(.19,-.17),(5.,8.),(.32,.09),(6.,9.),(.23,-.24),(8.,3.),(.18,.01), 00866 R (9.,4.) / 00867 C 00868 DATA CTRUE6/ 00869 A(.1,.1),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.),(1.,2.), 00870 B(.09,-.12),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.),(3.,4.), 00871 C (3.,4.), 00872 D(.03,-.09),(.15,-.03),(5.,6.),(5.,6.),(5.,6.),(5.,6.),(5.,6.), 00873 E (5.,6.), 00874 F(.03,.03),(-.18,.03),(.03,-.09),(7.,8.),(7.,8.),(7.,8.),(7.,8.), 00875 G (7.,8.), 00876 H(.09,.03),(.03,.12),(.12,.03),(.03,.06),(2.,3.),(2.,3.),(2.,3.), 00877 I (2.,3.), 00878 J(.1,.1),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.),(4.,5.), 00879 K(.09,-.12),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.),(6.,7.), 00880 L (6.,7.), 00881 M(.03,-.09),(8.,9.),(.15,-.03),(2.,5.),(2.,5.),(2.,5.),(2.,5.), 00882 N (2.,5.), 00883 O(.03,.03),(3.,6.),(-.18,.03),(4.,7.),(.03,-.09),(7.,2.),(7.,2.), 00884 P (7.,2.), 00885 Q(.09,.03),(5.,8.),(.03,.12),(6.,9.),(.12,.03),(8.,3.),(.03,.06), 00886 R (9.,4.) / 00887 C 00888 C 00889 DATA ITRUE2/ 0, 1, 2, 2, 3/ 00890 DATA ITRUE3/ 0, 1, 2, 2, 2/ 00891 C 00892 JUMP=ICASE-25 00893 DO 520 INCX=1,2 00894 DO 500 NP1=1,5 00895 N=NP1-1 00896 LEN= 2*MAX0(N,1) 00897 C SET VECTOR ARGUMENTS. 00898 DO 22 I=1,LEN 00899 SX(I) = SNGL(DV(I,NP1,INCX)) 00900 DX(I)=DV(I,NP1,INCX) 00901 22 CX(I)=CV(I,NP1,INCX) 00902 C 00903 C BRANCH TO INVOKE SUBPROGRAM TO BE TESTED. 00904 C 00905 GO TO (260,270,280,290,300,310,320, 00906 * 330,340,350,360,370,380),JUMP 00907 C 26. SNRM2 00908 260 STEMP=SNGL(DTRUE1(NP1)) 00909 CALL STEST1 (SNRM2(N,SX,INCX),STEMP,STEMP,SFAC) 00910 GO TO 500 00911 C 27. DNRM2 00912 270 CALL DTEST1 (DNRM2(N,DX,INCX),DTRUE1(NP1),DTRUE1(NP1),DFAC) 00913 GO TO 500 00914 C 28. SCNRM2 00915 280 CALL STEST1 (SCNRM2(N,CX,INCX),STRUE2(NP1),STRUE2(NP1),SFAC) 00916 GO TO 500 00917 C 29. SASUM 00918 290 STEMP=SNGL(DTRUE3(NP1)) 00919 CALL STEST1 (SASUM(N,SX,INCX),STEMP,STEMP,SFAC) 00920 GO TO 500 00921 C 30. DASUM 00922 300 CALL DTEST1 (DASUM(N,DX,INCX),DTRUE3(NP1),DTRUE3(NP1),DFAC) 00923 GO TO 500 00924 C 31. SCASUM 00925 310 CALL STEST1 (SCASUM(N,CX,INCX),STRUE4(NP1),STRUE4(NP1),SFAC) 00926 GO TO 500 00927 C 32. SSCALE 00928 320 CALL SSCAL(N,SA,SX,INCX) 00929 DO 322 I=1,LEN 00930 322 STRUE(I)= SNGL(DTRUE5(I,NP1,INCX)) 00931 CALL STEST(LEN,SX,STRUE,STRUE,SFAC) 00932 GO TO 500 00933 C 33. DSCALE 00934 330 CALL DSCAL(N,DA,DX,INCX) 00935 CALL DTEST(LEN,DX,DTRUE5(1,NP1,INCX),DTRUE5(1,NP1,INCX),DFAC) 00936 GO TO 500 00937 C 34. CSCALE 00938 340 CALL CSCAL(N,CA,CX,INCX) 00939 CALL CTEST (LEN,CX,CTRUE5(1,NP1,INCX),CTRUE5(1,NP1,INCX),SFAC) 00940 GO TO 500 00941 C 35. CSSCAL 00942 350 CALL CSSCAL(N,SA,CX,INCX) 00943 CALL CTEST (LEN,CX,CTRUE6(1,NP1,INCX),CTRUE6(1,NP1,INCX),SFAC) 00944 GO TO 500 00945 C 36. ISAMAX 00946 360 CALL ITEST1 (ISAMAX(N,SX,INCX),ITRUE2(NP1)) 00947 GO TO 500 00948 C 37. IDAMAX 00949 370 CALL ITEST1 (IDAMAX(N,DX,INCX),ITRUE2(NP1)) 00950 GO TO 500 00951 C 38. ICAMAX 00952 380 CALL ITEST1 (ICAMAX(N,CX,INCX),ITRUE3(NP1)) 00953 C 00954 500 CONTINUE 00955 520 CONTINUE 00956 RETURN 00957 END 00958 SUBROUTINE CHECK2(SFAC,SDFAC,DFAC,DQFAC) 00959 C1 ********************************* CHECK2 ************************* 00960 C THIS SUBPROGRAM TESTS THE BASIC LINEAR ALGEBRA SUBPROGRAMS 1-11, 00961 C 14-15, AND 18-25. SUBPROGRAMS IN THIS SET EACH REQUIRE TWO ARRAYS 00962 C IN THE PARAMETER LIST. 00963 C 00964 C C. L. LAWSON, JPL, 1975 FEB 26, APR 29, MAY 8, MAY 28 00965 C2 00966 COMMON/COMBLA/NPRINT,ICASE,N,INCX,INCY,MODE,PASS 00967 C 00968 LOGICAL PASS 00969 INTEGER INCXS(4),INCYS(4),LENS(4,2),NS(4),QC(10) 00970 REAL SX(7),SY(7),STX(7),STY(7),SSIZE1(4),SSIZE2(14,2) 00971 REAL SSIZE(7),SPARAM(5),ST7B(4,4),SSIZE3(4) 00972 DOUBLE PRECISION DX(7),DA,DX1(7),DY1(7),DY(7),DT7(4,4),DT8(7,4,4) 00973 DOUBLE PRECISION DX2(7), DY2(7), DT2(4,4,2), DPARAM(5), DPAR(5,4) 00974 DOUBLE PRECISION DSDOT,DDOT,DQDOTI,DQDOTA,DFAC,DQFAC 00975 DOUBLE PRECISION DT10X(7,4,4),DT10Y(7,4,4),DB 00976 DOUBLE PRECISION DSIZE1(4),DSIZE2(7,2),DSIZE(7) 00977 DOUBLE PRECISION DC,DS,DT9X(7,4,4),DT9Y(7,4,4),DTX(7),DTY(7) 00978 DOUBLE PRECISION DT19X(7,4,16),DT19XA(7,4,4),DT19XB(7,4,4) 00979 DOUBLE PRECISION DT19XC(7,4,4),DT19XD(7,4,4),DT19Y(7,4,16) 00980 DOUBLE PRECISION DT19YA(7,4,4),DT19YB(7,4,4),DT19YC(7,4,4) 00981 DOUBLE PRECISION DT19YD(7,4,4) 00982 C 00983 COMPLEX CX(7),CA,CX1(7),CY1(7),CY(7),CT6(4,4),CT7(4,4) 00984 COMPLEX CT8(7,4,4),CSIZE1(4),CSIZE2(7,2) 00985 COMPLEX CT10X(7,4,4), CT10Y(7,4,4) 00986 COMPLEX CDOT(1) 00987 COMPLEX CDOTC,CDOTU 00988 EQUIVALENCE (DT19X(1,1,1),DT19XA(1,1,1)),(DT19X(1,1,5), 00989 A DT19XB(1,1,1)),(DT19X(1,1,9),DT19XC(1,1,1)), 00990 B (DT19X(1,1,13),DT19XD(1,1,1)) 00991 EQUIVALENCE (DT19Y(1,1,1),DT19YA(1,1,1)),(DT19Y(1,1,5), 00992 A DT19YB(1,1,1)),(DT19Y(1,1,9),DT19YC(1,1,1)), 00993 B (DT19Y(1,1,13),DT19YD(1,1,1)) 00994 DATA SA,DA,CA,DB,SB/.3,.3D0,(.4,-.7),.25D0,.1/ 00995 DATA INCXS/ 1, 2, -2, -1 / 00996 DATA INCYS/ 1, -2, 1, -2 / 00997 DATA LENS/1, 1, 2, 4, 1, 1, 3, 7/ 00998 DATA NS / 0, 1, 2, 4 / 00999 DATA SC,SS,DC,DS/ .8,.6,.8D0,.6D0/ 01000 DATA DX1/ .6D0, .1D0,-.5D0, .8D0, .9D0,-.3D0,-.4D0/ 01001 DATA DY1/ .5D0,-.9D0, .3D0, .7D0,-.6D0, .2D0, .8D0/ 01002 DATA DX2/ 1.D0,.01D0, .02D0,1.25D0,.06D0, 2.D0, 1.D0/ 01003 DATA DY2/ 1.D0,.04D0,-.03D0,-1.D0,.05D0,3.D0,-1.D0/ 01004 DATA CX1/(.7,-.8),(-.4,-.7),(-.1,-.9),(.2,-.8),(-.9,-.4),(.1,.4), 01005 * (-.6,.6)/ 01006 DATA CY1/(.6,-.6),(-.9,.5),(.7,-.6),(.1,-.5),(-.1,-.2),(-.5,-.3), 01007 * (.8,-.7) / 01008 C 01009 C FOR DQDOTI AND DQDOTA 01010 C 01011 DATA DT2/0.25D0,1.25D0,1.2504D0,-0.0002D0, 01012 A 0.25D0,1.25D0,0.24D0,0.2492D0, 01013 B 0.25D0,1.25D0,0.31D0,0.2518D0, 01014 C 0.25D0,1.25D0,1.2497D0,0.0007D0, 01015 D 0.D0,2.D0,2.0008D0,-.5004D0, 01016 E 0.D0,2.D0,-.02D0,-.0016D0, 01017 F 0.D0,2.D0,.12D0,.0036D0, 01018 G 0.D0,2.D0,1.9994D0,-0.4986D0/ 01019 DATA DT7/ 0.D0,.30D0,.21D0,.62D0, 0.D0,.30D0,-.07D0,.85D0, 01020 * 0.D0,.30D0,-.79D0,-.74D0, 0.D0,.30D0,.33D0,1.27D0/ 01021 DATA ST7B/ .1, .4, .31, .72, .1, .4, .03, .95, 01022 * .1, .4, -.69, -.64, .1, .4, .43, 1.37/ 01023 C 01024 C FOR CDOTU 01025 C 01026 DATA CT7/(0.,0.),(-.06,-.90),(.65,-.47),(-.34,-1.22), 01027 1 (0.,0.),(-.06,-.90),(-.59,-1.46),(-1.04,-.04), 01028 2 (0.,0.),(-.06,-.90),(-.83,.59), ( .07,-.37), 01029 3 (0.,0.),(-.06,-.90),(-.76,-1.15),(-1.33,-1.82)/ 01030 C 01031 C FOR CDOTC 01032 C 01033 DATA CT6/(0.,0.),(.90,0.06), (.91,-.77), (1.80,-.10), 01034 A (0.,0.),(.90,0.06), (1.45,.74), (.20,.90), 01035 B (0.,0.),(.90,0.06), (-.55,.23), (.83,-.39), 01036 C (0.,0.),(.90,0.06), (1.04,0.79), (1.95,1.22)/ 01037 C 01038 DATA DT8/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01039 1 .68D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01040 2 .68D0,-.87D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01041 3 .68D0,-.87D0,.15D0,.94D0, 0.D0,0.D0,0.D0, 01042 4 .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01043 5 .68D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01044 6 .35D0,-.9D0,.48D0, 0.D0,0.D0,0.D0,0.D0, 01045 7 .38D0,-.9D0,.57D0,.7D0,-.75D0,.2D0,.98D0, 01046 8 .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01047 9 .68D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01048 A .35D0,-.72D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01049 B .38D0,-.63D0,.15D0,.88D0, 0.D0,0.D0,0.D0, 01050 C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01051 D .68D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01052 E .68D0,-.9D0,.33D0, 0.D0,0.D0,0.D0,0.D0, 01053 F .68D0,-.9D0,.33D0,.7D0,-.75D0,.2D0,1.04D0/ 01054 C 01055 DATA CT8/ 01056 A(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01057 B(.32,-1.41),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01058 C(.32,-1.41),(-1.55,.5),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01059 D(.32,-1.41),(-1.55,.5),(.03,-.89),(-.38,-.96),(0.,0.),(0.,0.), 01060 E (0.,0.), 01061 F(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01062 G(.32,-1.41),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01063 H(-.07,-.89),(-.9,.5),(.42,-1.41),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01064 I(.78,.06),(-.9,.5),(.06,-.13),(.1,-.5),(-.77,-.49),(-.5,-.3), 01065 J (.52,-1.51), 01066 K(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01067 L(.32,-1.41),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01068 M(-.07,-.89),(-1.18,-.31),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01069 N(.78,.06),(-1.54,.97),(.03,-.89),(-.18,-1.31),(0.,0.),(0.,0.), 01070 O(0.,0.),(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01071 P(.32,-1.41),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01072 Q(.32,-1.41),(-.9,.5),(.05,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01073 R(.32,-1.41),(-.9,.5),(.05,-.6),(.1,-.5),(-.77,-.49),(-.5,-.3), 01074 S (.32,-1.16) / 01075 C 01076 C 01077 C TRUE X VALUES AFTER ROTATION USING SROT OR DROT. 01078 DATA DT9X/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01079 A .78D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01080 B .78D0,-.46D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01081 C .78D0,-.46D0,-.22D0,1.06D0, 0.D0,0.D0,0.D0, 01082 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01083 E .78D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01084 F .66D0,.1D0,-.1D0, 0.D0,0.D0,0.D0,0.D0, 01085 G .96D0,.1D0,-.76D0,.8D0,.90D0,-.3D0,-.02D0, 01086 H .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01087 I .78D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01088 J -.06D0,.1D0,-.1D0, 0.D0,0.D0,0.D0,0.D0, 01089 K .90D0,.1D0,-.22D0,.8D0,.18D0,-.3D0,-.02D0, 01090 L .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01091 M .78D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01092 N .78D0,.26D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01093 O .78D0,.26D0,-.76D0,1.12D0, 0.D0,0.D0,0.D0/ 01094 C 01095 C TRUE Y VALUES AFTER ROTATION USING SROT OR DROT. 01096 C 01097 DATA DT9Y/ .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01098 A .04D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01099 B .04D0,-.78D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01100 C .04D0,-.78D0, .54D0, .08D0, 0.D0,0.D0,0.D0, 01101 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01102 E .04D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01103 F .7D0,-.9D0,-.12D0, 0.D0,0.D0,0.D0,0.D0, 01104 G .64D0,-.9D0,-.30D0, .7D0,-.18D0, .2D0, .28D0, 01105 H .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01106 I .04D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01107 J .7D0,-1.08D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01108 K .64D0,-1.26D0,.54D0, .20D0, 0.D0,0.D0,0.D0, 01109 L .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01110 M .04D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01111 N .04D0,-.9D0, .18D0, 0.D0,0.D0,0.D0,0.D0, 01112 O .04D0,-.9D0, .18D0, .7D0,-.18D0, .2D0, .16D0/ 01113 C 01114 DATA DT10X/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01115 A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01116 B .5D0,-.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01117 C .5D0,-.9D0,.3D0,.7D0, 0.D0,0.D0,0.D0, 01118 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01119 E .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01120 F .3D0,.1D0 ,.5D0, 0.D0,0.D0,0.D0,0.D0, 01121 G .8D0,.1D0 ,-.6D0,.8D0 ,.3D0,-.3D0,.5D0, 01122 H .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01123 I .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01124 J -.9D0,.1D0,.5D0, 0.D0,0.D0,0.D0,0.D0, 01125 K .7D0, .1D0,.3D0, .8D0,-.9D0,-.3D0,.5D0, 01126 L .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01127 M .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01128 N .5D0,.3D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01129 O .5D0,.3D0,-.6D0,.8D0, 0.D0,0.D0,0.D0/ 01130 C 01131 DATA DT10Y/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01132 A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01133 B .6D0,.1D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01134 C .6D0,.1D0,-.5D0,.8D0, 0.D0,0.D0,0.D0, 01135 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01136 E .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01137 F -.5D0,-.9D0,.6D0, 0.D0,0.D0,0.D0,0.D0, 01138 G -.4D0,-.9D0,.9D0, .7D0,-.5D0, .2D0,.6D0, 01139 H .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01140 I .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01141 J -.5D0,.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01142 K -.4D0,.9D0,-.5D0,.6D0, 0.D0,0.D0,0.D0, 01143 L .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01144 M .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01145 N .6D0,-.9D0,.1D0, 0.D0,0.D0,0.D0,0.D0, 01146 O .6D0,-.9D0,.1D0, .7D0,-.5D0, .2D0, .8D0/ 01147 C 01148 DATA CT10X/ 01149 A(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01150 B(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01151 C(.6,-.6),(-.9,.5),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01152 D(.6,-.6),(-.9,.5),(.7,-.6),(.1,-.5),(0.,0.),(0.,0.),(0.,0.), 01153 E(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01154 F(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01155 G(.7,-.6),(-.4,-.7),(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01156 H(.8,-.7),(-.4,-.7),(-.1,-.2),(.2,-.8),(.7,-.6),(.1,.4),(.6,-.6), 01157 I(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01158 J(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01159 K(-.9,.5),(-.4,-.7),(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01160 L(.1,-.5),(-.4,-.7),(.7,-.6),(.2,-.8),(-.9,.5),(.1,.4),(.6,-.6), 01161 M(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01162 N(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01163 O(.6,-.6),(.7,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01164 P(.6,-.6),(.7,-.6),(-.1,-.2),(.8,-.7),(0.,0.),(0.,0.),(0.,0.) / 01165 C 01166 DATA CT10Y/ 01167 A(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01168 B(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01169 C(.7,-.8),(-.4,-.7),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01170 D(.7,-.8),(-.4,-.7),(-.1,-.9),(.2,-.8),(0.,0.),(0.,0.),(0.,0.), 01171 E(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01172 F(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01173 G(-.1,-.9),(-.9,.5),(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01174 H(-.6,.6),(-.9,.5),(-.9,-.4),(.1,-.5),(-.1,-.9),(-.5,-.3),(.7,-.8), 01175 I(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01176 J(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01177 K(-.1,-.9),(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01178 L(-.6,.6),(-.9,-.4),(-.1,-.9),(.7,-.8),(0.,0.),(0.,0.),(0.,0.), 01179 M(.6,-.6),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01180 N(.7,-.8),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01181 O(.7,-.8),(-.9,.5),(-.4,-.7),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01182 P(.7,-.8),(-.9,.5),(-.4,-.7),(.1,-.5),(-.1,-.9),(-.5,-.3),(.2,-.8)/ 01183 C TRUE X RESULTS F0R ROTATIONS SROTM AND DROTM 01184 DATA DT19XA/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01185 A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01186 B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01187 C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01188 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01189 E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01190 F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01191 G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01192 H .6D0, .1D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01193 I -.8D0, 3.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01194 J -.9D0, 2.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01195 K 3.5D0, -.4D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01196 L .6D0, .1D0, -.5D0, .8D0, 0.D0,0.D0,0.D0, 01197 M -.8D0, 3.8D0, -2.2D0, -1.2D0, 0.D0,0.D0,0.D0, 01198 N -.9D0, 2.8D0, -1.4D0, -1.3D0, 0.D0,0.D0,0.D0, 01199 O 3.5D0, -.4D0, -2.2D0, 4.7D0, 0.D0,0.D0,0.D0/ 01200 C 01201 DATA DT19XB/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01202 A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01203 B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01204 C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01205 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01206 E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01207 F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01208 G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01209 H .6D0, .1D0, -.5D0, 0.D0,0.D0,0.D0,0.D0, 01210 I 0.D0, .1D0, -3.0D0, 0.D0,0.D0,0.D0,0.D0, 01211 J -.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, 01212 K 3.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, 01213 L .6D0, .1D0, -.5D0, .8D0, .9D0, -.3D0, -.4D0, 01214 M -2.0D0, .1D0, 1.4D0, .8D0, .6D0, -.3D0, -2.8D0, 01215 N -1.8D0, .1D0, 1.3D0, .8D0, 0.D0, -.3D0, -1.9D0, 01216 O 3.8D0, .1D0, -3.1D0, .8D0, 4.8D0, -.3D0, -1.5D0 / 01217 C 01218 DATA DT19XC/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01219 A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01220 B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01221 C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01222 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01223 E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01224 F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01225 G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01226 H .6D0, .1D0, -.5D0, 0.D0,0.D0,0.D0,0.D0, 01227 I 4.8D0, .1D0, -3.0D0, 0.D0,0.D0,0.D0,0.D0, 01228 J 3.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, 01229 K 2.1D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, 01230 L .6D0, .1D0, -.5D0, .8D0, .9D0, -.3D0, -.4D0, 01231 M -1.6D0, .1D0, -2.2D0, .8D0, 5.4D0, -.3D0, -2.8D0, 01232 N -1.5D0, .1D0, -1.4D0, .8D0, 3.6D0, -.3D0, -1.9D0, 01233 O 3.7D0, .1D0, -2.2D0, .8D0, 3.6D0, -.3D0, -1.5D0 / 01234 C 01235 DATA DT19XD/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01236 A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01237 B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01238 C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01239 D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01240 E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01241 F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01242 G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01243 H .6D0, .1D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01244 I -.8D0, -1.0D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01245 J -.9D0, -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01246 K 3.5D0, .8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01247 L .6D0, .1D0, -.5D0, .8D0, 0.D0,0.D0,0.D0, 01248 M -.8D0, -1.0D0, 1.4D0, -1.6D0, 0.D0,0.D0,0.D0, 01249 N -.9D0, -.8D0, 1.3D0, -1.6D0, 0.D0,0.D0,0.D0, 01250 O 3.5D0, .8D0, -3.1D0, 4.8D0, 0.D0,0.D0,0.D0/ 01251 C TRUE Y RESULTS FOR ROTATIONS SROTM AND DROTM 01252 DATA DT19YA/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01253 A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01254 B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01255 C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01256 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01257 E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01258 F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01259 G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01260 H .5D0, -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01261 I .7D0, -4.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01262 J 1.7D0, -.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01263 K -2.6D0, 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01264 L .5D0, -.9D0, .3D0, .7D0, 0.D0,0.D0,0.D0, 01265 M .7D0, -4.8D0, 3.0D0, 1.1D0, 0.D0,0.D0,0.D0, 01266 N 1.7D0, -.7D0, -.7D0, 2.3D0, 0.D0,0.D0,0.D0, 01267 O -2.6D0, 3.5D0, -.7D0, -3.6D0, 0.D0,0.D0,0.D0/ 01268 C 01269 DATA DT19YB/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01270 A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01271 B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01272 C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01273 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01274 E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01275 F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01276 G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01277 H .5D0, -.9D0, .3D0, 0.D0,0.D0,0.D0,0.D0, 01278 I 4.0D0, -.9D0, -.3D0, 0.D0,0.D0,0.D0,0.D0, 01279 J -.5D0, -.9D0, 1.5D0, 0.D0,0.D0,0.D0,0.D0, 01280 K -1.5D0, -.9D0, -1.8D0, 0.D0,0.D0,0.D0,0.D0, 01281 L .5D0, -.9D0, .3D0, .7D0, -.6D0, .2D0, .8D0, 01282 M 3.7D0, -.9D0, -1.2D0, .7D0, -1.5D0, .2D0, 2.2D0, 01283 N -.3D0, -.9D0, 2.1D0, .7D0, -1.6D0, .2D0, 2.0D0, 01284 O -1.6D0, -.9D0, -2.1D0, .7D0, 2.9D0, .2D0, -3.8D0 / 01285 C 01286 DATA DT19YC/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01287 A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01288 B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01289 C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01290 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01291 E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01292 F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01293 G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01294 H .5D0, -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01295 I 4.0D0, -6.3D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01296 J -.5D0, .3D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01297 K -1.5D0, 3.0D0, 0.D0,0.D0,0.D0,0.D0,0.D0, 01298 L .5D0, -.9D0, .3D0, .7D0, 0.D0,0.D0,0.D0, 01299 M 3.7D0, -7.2D0, 3.0D0, 1.7D0, 0.D0,0.D0,0.D0, 01300 N -.3D0, .9D0, -.7D0, 1.9D0, 0.D0,0.D0,0.D0, 01301 O -1.6D0, 2.7D0, -.7D0, -3.4D0, 0.D0,0.D0,0.D0/ 01302 C 01303 DATA DT19YD/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01304 A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01305 B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01306 C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01307 D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01308 E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01309 F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01310 G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01311 H .5D0, -.9D0, .3D0, 0.D0,0.D0,0.D0,0.D0, 01312 I .7D0, -.9D0, 1.2D0, 0.D0,0.D0,0.D0,0.D0, 01313 J 1.7D0, -.9D0, .5D0, 0.D0,0.D0,0.D0,0.D0, 01314 K -2.6D0, -.9D0, -1.3D0, 0.D0,0.D0,0.D0,0.D0, 01315 L .5D0, -.9D0, .3D0, .7D0, -.6D0, .2D0, .8D0, 01316 M .7D0, -.9D0, 1.2D0, .7D0, -1.5D0, .2D0, 1.6D0, 01317 N 1.7D0, -.9D0, .5D0, .7D0, -1.6D0, .2D0, 2.4D0, 01318 O -2.6D0, -.9D0, -1.3D0, .7D0, 2.9D0, .2D0, -4.0D0 / 01319 C 01320 DATA SSIZE1/ 0. , .3 , 1.6 , 3.2 / 01321 DATA DSIZE1/ 0.D0, .3D0, 1.6D0, 3.2D0 / 01322 DATA SSIZE3/ .1, .4, 1.7, 3.3 / 01323 C 01324 C FOR CDOTC AND CDOTU 01325 C 01326 DATA CSIZE1/ (0.,0.), (.9,.9), (1.63,1.73), (2.90,2.78) / 01327 DATA SSIZE2/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 01328 A 1.17,1.17,1.17,1.17,1.17,1.17,1.17, 01329 B 1.17,1.17,1.17,1.17,1.17,1.17,1.17/ 01330 DATA DSIZE2/0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 01331 A 1.17D0,1.17D0,1.17D0,1.17D0,1.17D0,1.17D0,1.17D0/ 01332 C 01333 C FOR CAXPY 01334 C 01335 DATA CSIZE2/ 01336 A (0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.), 01337 B (1.54,1.54),(1.54,1.54),(1.54,1.54),(1.54,1.54),(1.54,1.54), 01338 C (1.54,1.54),(1.54,1.54) / 01339 C 01340 C FOR SROTM AND DROTM 01341 C 01342 DATA DPAR/-2.D0, 0.D0,0.D0,0.D0,0.D0, 01343 A -1.D0, 2.D0, -3.D0, -4.D0, 5.D0, 01344 B 0.D0, 0.D0, 2.D0, -3.D0, 0.D0, 01345 C 1.D0, 5.D0, 2.D0, 0.D0, -4.D0/ 01346 C 01347 DO 520 KI=1,4 01348 INCX = INCXS(KI) 01349 INCY = INCYS(KI) 01350 MX = IABS(INCX) 01351 MY = IABS(INCY) 01352 C 01353 DO 500 KN=1,4 01354 N= NS(KN) 01355 KSIZE=MIN0(2,KN) 01356 LENX = LENS(KN,MX) 01357 LENY = LENS(KN,MY) 01358 C INITIALIZE ALL ARGUMENT ARRAYS. 01359 DO 5 I=1,7 01360 SX(I)= SNGL(DX1(I)) 01361 SY(I)= SNGL(DY1(I)) 01362 DX(I)= DX1(I) 01363 DY(I)= DY1(I) 01364 CX(I)= CX1(I) 01365 5 CY(I)= CY1(I) 01366 C 01367 C BRANCH TO SELECT SUBPROGRAM TO BE TESTED. 01368 C 01369 GO TO ( 10, 20, 30, 40, 50, 60, 70, 80, 90,100, 01370 A 110,999,999,140,150,999,999,180,190,200, 01371 B 210,220,230,240,250), ICASE 01372 C 1. SDOT 01373 10 CALL STEST1 (SDOT(N,SX,INCX,SY,INCY),SNGL(DT7(KN,KI)), 01374 * SSIZE1(KN),SFAC) 01375 GO TO 500 01376 C 2. DSDOT 01377 20 CALL STEST1 (SNGL(DSDOT(N,SX,INCX,SY,INCY)), 01378 * SNGL(DT7(KN,KI)),SSIZE1(KN),SFAC) 01379 GO TO 500 01380 C 3. SDSDOT 01381 30 CALL STEST1 (SDSDOT(N,SB,SX,INCX,SY,INCY), 01382 * ST7B(KN,KI),SSIZE3(KN),SFAC) 01383 GO TO 500 01384 C 4. DDOT 01385 40 CALL DTEST1 (DDOT(N,DX,INCX,DY,INCY),DT7(KN,KI), 01386 * DSIZE1(KN),DFAC) 01387 GO TO 500 01388 C 5. DQDOTI 01389 50 CONTINUE 01390 C DQDOTI AND DQDOTA ARE SUPPOSED TO USE EXTENDED 01391 C PRECISION ARITHMETIC INTERNALLY. 01392 C SET MODE = 1 OR 2 TO DISTINGUISH TESTS OF DQDOTI OR DQDOTA 01393 C IN THE DIAGNOSTIC OUTPUT. 01394 C 01395 MODE = 1 01396 CALL DTEST1 (DQDOTI(N,DB,QC,DX2,INCX,DY2,INCY), 01397 * DT2(KN,KI,1),DT2(KN,KI,1),DQFAC) 01398 GO TO 500 01399 C 6. DQDOTA 01400 60 CONTINUE 01401 C TO TEST DQDOTA WE ACTUALLY TEST BOTH DQDOTI AND DQDOTA. 01402 C THE OUTPUT VALUE OF QX FROM DQDOTI WILL BE USED AS INPUT 01403 C TO DQDOTA. QX IS SUPPOSED TO BE IN A MACHINE-DEPENDENT 01404 C EXTENDED PRECISION FORM. 01405 C MODE IS SET TO 1 OR 2 TO DISTINGUISH TESTS OF 01406 C DQDOTI OR DQDOTA IN THE DIAGNOSTIC OUTPUT. 01407 C 01408 MODE = 1 01409 CALL DTEST1 (DQDOTI(N,DB,QC,DX2,INCX,DY2,INCY), 01410 * DT2(KN,KI,1),DT2(KN,KI,1),DQFAC) 01411 MODE = 2 01412 CALL DTEST1 (DQDOTA(N,-DB,QC,DX2,INCX,DY2,INCY), 01413 * DT2(KN,KI,2),DT2(KN,KI,2),DQFAC) 01414 GO TO 500 01415 C 7. CDOTC 01416 70 CDOT(1) = CDOTC(N,CX,INCX,CY,INCY) 01417 CALL CTEST(1,CDOT,CT6(KN,KI),CSIZE1(KN),SFAC) 01418 GO TO 500 01419 C 8. CDOTU 01420 80 CDOT(1) = CDOTU(N,CX,INCX,CY,INCY) 01421 CALL CTEST(1,CDOT,CT7(KN,KI),CSIZE1(KN),SFAC) 01422 GO TO 500 01423 C 9. SAXPY 01424 90 CALL SAXPY(N,SA,SX,INCX,SY,INCY) 01425 DO 95 J=1,LENY 01426 95 STY(J)= SNGL(DT8(J,KN,KI)) 01427 CALL STEST(LENY,SY,STY,SSIZE2(1,KSIZE),SFAC) 01428 GO TO 500 01429 C 10. DAXPY 01430 100 CALL DAXPY(N,DA,DX,INCX,DY,INCY) 01431 CALL DTEST(LENY,DY,DT8(1,KN,KI),DSIZE2(1,KSIZE),DFAC) 01432 GO TO 500 01433 C 11. CAXPY 01434 110 CALL CAXPY(N,CA,CX,INCX,CY,INCY) 01435 CALL CTEST (LENY,CY,CT8(1,KN,KI),CSIZE2(1,KSIZE),SFAC) 01436 GO TO 500 01437 C 14. SROT 01438 140 CONTINUE 01439 DO 144 I=1,7 01440 SX(I)= SNGL(DX1(I)) 01441 SY(I)= SNGL(DY1(I)) 01442 STX(I)= SNGL(DT9X(I,KN,KI)) 01443 STY(I)= SNGL(DT9Y(I,KN,KI)) 01444 144 CONTINUE 01445 CALL SROT (N,SX,INCX,SY,INCY,SC,SS) 01446 CALL STEST(LENX,SX,STX,SSIZE2(1,KSIZE),SFAC) 01447 CALL STEST(LENY,SY,STY,SSIZE2(1,KSIZE),SFAC) 01448 GO TO 500 01449 C 15. DROT 01450 150 CONTINUE 01451 DO 154 I=1,7 01452 DX(I)=DX1(I) 01453 DY(I)=DY1(I) 01454 154 CONTINUE 01455 CALL DROT (N,DX,INCX,DY,INCY,DC,DS) 01456 CALL DTEST(LENX,DX,DT9X(1,KN,KI), DSIZE2(1,KSIZE),DFAC) 01457 CALL DTEST(LENY,DY,DT9Y(1,KN,KI),DSIZE2(1,KSIZE),DFAC) 01458 GO TO 500 01459 C 18. SROTM 01460 180 KNI=KN+4*(KI-1) 01461 DO 189 KPAR=1,4 01462 DO 182 I=1,7 01463 SX(I)= SNGL(DX1(I)) 01464 SY(I)= SNGL(DY1(I)) 01465 STX(I)= SNGL(DT19X(I,KPAR,KNI)) 01466 182 STY(I)= SNGL(DT19Y(I,KPAR,KNI)) 01467 C 01468 DO 186 I=1,5 01469 186 SPARAM(I) = SNGL(DPAR(I,KPAR)) 01470 C SET MODE TO IDENTIFY DIAGNOSTIC OUTPUT, 01471 C IF ANY 01472 MODE = INT(SPARAM(1)) 01473 C 01474 DO 187 I=1,LENX 01475 187 SSIZE(I)=STX(I) 01476 C THE TRUE RESULTS DT19X(1,2,7) AND 01477 C DT19X(5,3,8) ARE ZERO DUE TO CANCELLATION. 01478 C DT19X(1,2,7) = 2.*.6 - 4.*.3 = 0 01479 C DT19X(5,3,8) = .9 - 3.*.3 = 0 01480 C FOR THESE CASES RESPECTIVELY SET SIZE( ) 01481 C EQUAL TO 2.4 AND 1.8 01482 IF ((KPAR .EQ. 2) .AND. (KNI .EQ. 7)) 01483 1 SSIZE(1) = 2.4E0 01484 IF ((KPAR .EQ. 3) .AND. (KNI .EQ. 8)) 01485 1 SSIZE(5) = 1.8E0 01486 C 01487 CALL SROTM(N,SX,INCX,SY,INCY,SPARAM) 01488 CALL STEST(LENX,SX,STX,SSIZE,SFAC) 01489 CALL STEST(LENY,SY,STY,STY,SFAC) 01490 189 CONTINUE 01491 GO TO 500 01492 C 19. DROTM 01493 190 KNI=KN+4*(KI-1) 01494 DO 199 KPAR=1,4 01495 DO 192 I=1,7 01496 DX(I)=DX1(I) 01497 DY(I)=DY1(I) 01498 DTX(I)= DT19X(I,KPAR,KNI) 01499 192 DTY(I)= DT19Y(I,KPAR,KNI) 01500 C 01501 DO 196 I=1,5 01502 196 DPARAM(I) = DPAR(I,KPAR) 01503 C SET MODE TO IDENTIFY DIAGNOSTIC OUTPUT, 01504 C IF ANY 01505 MODE = IDINT(DPARAM(1)) 01506 C 01507 DO 197 I=1,LENX 01508 197 DSIZE(I)=DTX(I) 01509 C SEE REMARK ABOVE ABOUT DT11X(1,2,7) 01510 C AND DT11X(5,3,8). 01511 IF ((KPAR .EQ. 2) .AND. (KNI .EQ. 7)) 01512 1 DSIZE(1) = 2.4D0 01513 IF ((KPAR .EQ. 3) .AND. (KNI .EQ. 8)) 01514 1 DSIZE(5) = 1.8D0 01515 C 01516 CALL DROTM(N,DX,INCX,DY,INCY,DPARAM) 01517 CALL DTEST(LENX,DX,DTX,DSIZE,DFAC) 01518 CALL DTEST(LENY,DY,DTY,DTY,DFAC) 01519 199 CONTINUE 01520 GO TO 500 01521 C 20. SCOPY 01522 200 DO 205 I=1,7 01523 205 STY(I)= SNGL(DT10Y( I,KN,KI)) 01524 CALL SCOPY(N,SX,INCX,SY,INCY) 01525 CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.) 01526 GO TO 500 01527 C 21. DCOPY 01528 210 CALL DCOPY(N,DX,INCX,DY,INCY) 01529 CALL DTEST(LENY,DY,DT10Y(1,KN,KI), DSIZE2(1,1),1.D0 ) 01530 GO TO 500 01531 C 22. CCOPY 01532 220 CALL CCOPY(N,CX,INCX,CY,INCY) 01533 CALL CTEST (LENY,CY,CT10Y(1,KN,KI),SSIZE2(1,1),1.) 01534 GO TO 500 01535 C 23. SSWAP 01536 230 CALL SSWAP(N,SX,INCX,SY,INCY) 01537 DO 235 I=1,7 01538 STX(I)= SNGL(DT10X(I,KN,KI)) 01539 235 STY(I)= SNGL(DT10Y(I,KN,KI)) 01540 CALL STEST(LENX,SX,STX,SSIZE2(1,1),1.) 01541 CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.) 01542 GO TO 500 01543 C 24. DSWAP 01544 240 CALL DSWAP(N,DX,INCX,DY,INCY) 01545 CALL DTEST(LENX,DX,DT10X(1,KN,KI), DSIZE2(1,1),1.D0) 01546 CALL DTEST(LENY,DY,DT10Y(1,KN,KI), DSIZE2(1,1),1.D0) 01547 GO TO 500 01548 C 25. CSWAP 01549 250 CALL CSWAP(N,CX,INCX,CY,INCY) 01550 CALL CTEST (LENX,CX,CT10X(1,KN,KI), SSIZE2(1,1),1.) 01551 CALL CTEST (LENY,CY, CT10Y(1,KN,KI), SSIZE2(1,1),1.) 01552 C 01553 C 01554 C 01555 500 CONTINUE 01556 520 CONTINUE 01557 RETURN 01558 C THE FOLLOWING STOP SHOULD NEVER BE REACHED. 01559 999 STOP 01560 END 01561 SUBROUTINE MPBLAS(I1) 01562 C 01563 C THIS SUBROUTINE IS CALLED TO SET UP BRENT'S MP PACKAGE 01564 C FOR USE BY THE EXTENDED PRECISION INNER PRODUCTS FROM THE BLAS. 01565 C 01566 C THE COMMON BLOCK FOR THE MP PACKAGE (MODIFIED TO GIVE IT A NAME) 01567 COMMON /MPCOM/ MPB, MPT, MPM, MPLUN, MPMXR, MPR(12) 01568 C 01569 C SET I1 = 1 TO FLAG THAT THIS ROUTINE WAS CALLED 01570 I1 = 1 01571 C 01572 C FOR FULL EXTENDED PRECISION ACCURACY, MPB SHOULD BE AS LARGE AS 01573 C POSSIBLE, SUBJECT TO THE RESTRICTIONS IN BRENT'S PAPER (ACM TRANS. 01574 C ON MATH. SOFTWARE, MARCH 1978, VOL. 4, NO. 1.). 01575 C STATEMENTS BELOW ARE FOR AN INTEGER WORDLENGTH OF 48, 36, 32, 01576 C 24, 18, AND 16. PICK ONE, OR GENERATE A NEW ONE. 01577 C 48 MPB = 4194304 01578 C 36 MPB = 65536 01579 32 MPB = 16384 01580 C 24 MPB = 1024 01581 C 18 MPB = 128 01582 C 16 MPB = 64 01583 C 01584 C SET UP REMAINING PARAMETERS 01585 C UNIT FOR ERROR MESSAGES 01586 MPLUN = 6 01587 C NUMBER OF MP DIGITS 01588 MPT = 8 01589 C DIMENSION OF R 01590 MPMXR = 12 01591 C EXPONENT RANGE 01592 MPM = 32767 01593 RETURN 01594 END 01595 C $$ ****** MPADD ****** 01596 SUBROUTINE MPADD (X, Y, Z) 01597 C ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP 01598 C NUMBERS. FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING. 01599 INTEGER X(1), Y(1), Z(1) 01600 CALL MPADD2 (X, Y, Z, Y, 0) 01601 RETURN 01602 END 01603 C $$ ****** MPMLP ****** 01604 SUBROUTINE MPMLP (U, V, W, J) 01605 C PERFORMS INNER MULTIPLICATION LOOP FOR MPMUL 01606 C NOTE THAT CARRIES ARE NOT PROPAGATED IN INNER LOOP, 01607 C WHICH SAVES TIME AT THE EXPENSE OF SPACE. 01608 INTEGER U(1), V(1), W 01609 DO 10 I = 1, J 01610 10 U(I) = U(I) + W*V(I) 01611 RETURN 01612 END 01613 C $$ ****** MPUNFL ****** 01614 SUBROUTINE MPUNFL (X) 01615 C CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE 01616 C EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M. 01617 INTEGER X(1) 01618 C SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC. 01619 CALL MPCHK (1, 4) 01620 C THE UNDERFLOWING NUMBER IS SET TO ZERO 01621 C AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN, 01622 C POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION 01623 C AFTER A PRESET NUMBER OF UNDERFLOWS. ACTION COULD EASILY 01624 C BE DETERMINED BY A FLAG IN LABELLED COMMON. 01625 X(1) = 0 01626 RETURN 01627 END 01628 C $$ ****** MPMULI ****** 01629 SUBROUTINE MPMULI (X, IY, Z) 01630 C MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z. 01631 C THIS IS FASTER THAN USING MPMUL. RESULT IS ROUNDED. 01632 C MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER 01633 C EVEN IF THE LAST DIGIT IS B. 01634 INTEGER X(1), Z(1) 01635 CALL MPMUL2 (X, IY, Z, 0) 01636 RETURN 01637 END 01638 C $$ ****** MPADD2 ****** 01639 SUBROUTINE MPADD2 (X, Y, Z, Y1, TRUNC) 01640 C MODIFIED FOR USE WITH BLAS. 01641 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01642 C CALLED BY MPADD, MPSUB ETC. 01643 C X, Y AND Z ARE MP NUMBERS, Y1 AND TRUNC ARE INTEGERS. 01644 C TO FORCE CALL BY REFERENCE RATHER THAN VALUE/RESULT, Y1 IS 01645 C DECLARED AS AN ARRAY, BUT ONLY Y1(1) IS EVER USED. 01646 C SETS Z = X + Y1(1)*ABS(Y), WHERE Y1(1) = +- Y(1). 01647 C IF TRUNC.EQ.0 R*-ROUNDING IS USED, OTHERWISE TRUNCATION. 01648 C R*-ROUNDING IS DEFINED IN KUKI AND CODI, COMM. ACM 01649 C 16(1973), 223. (SEE ALSO BRENT, IEEE TC-22(1973), 601.) 01650 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01651 INTEGER B, T, R, X(1), Y(1), Z(1), Y1(1), TRUNC 01652 INTEGER S, ED, RS, RE 01653 C CHECK FOR X OR Y ZERO 01654 IF (X(1).NE.0) GO TO 20 01655 10 CALL MPSTR(Y, Z) 01656 Z(1) = Y1(1) 01657 RETURN 01658 20 IF (Y1(1).NE.0) GO TO 40 01659 30 CALL MPSTR (X, Z) 01660 RETURN 01661 C COMPARE SIGNS 01662 40 S = X(1)*Y1(1) 01663 IF (IABS(S).LE.1) GO TO 60 01664 CALL MPCHK (1, 4) 01665 WRITE (LUN, 50) 01666 50 FORMAT (44H *** SIGN NOT 0, +1 OR -1 IN CALL TO MPADD2,, 01667 $ 33H POSSIBLE OVERWRITING PROBLEM ***) 01668 CALL MPERR 01669 Z(1) = 0 01670 RETURN 01671 C COMPARE EXPONENTS 01672 60 ED = X(2) - Y(2) 01673 MED = IABS(ED) 01674 IF (ED) 90, 70, 120 01675 C EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. 01676 70 IF (S.GT.0) GO TO 100 01677 DO 80 J = 1, T 01678 IF (X(J+2) - Y(J+2)) 100, 80, 130 01679 80 CONTINUE 01680 C RESULT IS ZERO 01681 Z(1) = 0 01682 RETURN 01683 C HERE EXPONENT(Y) .GE. EXPONENT(X) 01684 90 IF (MED.GT.T) GO TO 10 01685 100 RS = Y1(1) 01686 RE = Y(2) 01687 CALL MPADD3 (X, Y, S, MED, RE) 01688 C NORMALIZE, ROUND OR TRUNCATE, AND RETURN 01689 110 CALL MPNZR (RS, RE, Z, TRUNC) 01690 RETURN 01691 C ABS(X) .GT. ABS(Y) 01692 120 IF (MED.GT.T) GO TO 30 01693 130 RS = X(1) 01694 RE = X(2) 01695 CALL MPADD3 (Y, X, S, MED, RE) 01696 GO TO 110 01697 END 01698 C $$ ****** MPADD3 ****** 01699 SUBROUTINE MPADD3 (X, Y, S, MED, RE) 01700 C MODIFIED FOR USE WITH BLAS. 01701 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01702 C CALLED BY MPADD2, DOES INNER LOOPS OF ADDITION 01703 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01704 INTEGER B, T, R, X(1), Y(1), S, RE, C, TED 01705 TED = T + MED 01706 I2 = T + 4 01707 I = I2 01708 C = 0 01709 C CLEAR GUARD DIGITS TO RIGHT OF X DIGITS 01710 10 IF (I.LE.TED) GO TO 20 01711 R(I) = 0 01712 I = I - 1 01713 GO TO 10 01714 20 IF (S.LT.0) GO TO 130 01715 C HERE DO ADDITION, EXPONENT(Y) .GE. EXPONENT(X) 01716 IF (I.LT.T) GO TO 40 01717 30 J = I - MED 01718 R(I) = X(J+2) 01719 I = I - 1 01720 IF (I.GT.T) GO TO 30 01721 40 IF (I.LE.MED) GO TO 60 01722 J = I - MED 01723 C = Y(I+2) + X(J+2) + C 01724 IF (C.LT.B) GO TO 50 01725 C CARRY GENERATED HERE 01726 R(I) = C - B 01727 C = 1 01728 I = I - 1 01729 GO TO 40 01730 C NO CARRY GENERATED HERE 01731 50 R(I) = C 01732 C = 0 01733 I = I - 1 01734 GO TO 40 01735 60 IF (I.LE.0) GO TO 90 01736 C = Y(I+2) + C 01737 IF (C.LT.B) GO TO 70 01738 R(I) = 0 01739 C = 1 01740 I = I - 1 01741 GO TO 60 01742 70 R(I) = C 01743 I = I - 1 01744 C NO CARRY POSSIBLE HERE 01745 80 IF (I.LE.0) RETURN 01746 R(I) = Y(I+2) 01747 I = I - 1 01748 GO TO 80 01749 90 IF (C.EQ.0) RETURN 01750 C MUST SHIFT RIGHT HERE AS CARRY OFF END 01751 I2P = I2 + 1 01752 DO 100 J = 2, I2 01753 I = I2P - J 01754 100 R(I+1) = R(I) 01755 R(1) = 1 01756 RE = RE + 1 01757 RETURN 01758 C HERE DO SUBTRACTION, ABS(Y) .GT. ABS(X) 01759 110 J = I - MED 01760 R(I) = C - X(J+2) 01761 C = 0 01762 IF (R(I).GE.0) GO TO 120 01763 C BORROW GENERATED HERE 01764 C = -1 01765 R(I) = R(I) + B 01766 120 I = I - 1 01767 130 IF (I.GT.T) GO TO 110 01768 140 IF (I.LE.MED) GO TO 160 01769 J = I - MED 01770 C = Y(I+2) + C - X(J+2) 01771 IF (C.GE.0) GO TO 150 01772 C BORROW GENERATED HERE 01773 R(I) = C + B 01774 C = -1 01775 I = I - 1 01776 GO TO 140 01777 C NO BORROW GENERATED HERE 01778 150 R(I) = C 01779 C = 0 01780 I = I - 1 01781 GO TO 140 01782 160 IF (I.LE.0) RETURN 01783 C = Y(I+2) + C 01784 IF (C.GE.0) GO TO 70 01785 R(I) = C + B 01786 C = -1 01787 I = I - 1 01788 GO TO 160 01789 END 01790 C $$ ****** MPCDM ****** 01791 SUBROUTINE MPCDM (DX, Z) 01792 C MODIFIED FOR USE WITH BLAS. 01793 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01794 C CONVERTS DOUBLE-PRECISION NUMBER DX TO MULTIPLE-PRECISION Z. 01795 C SOME NUMBERS WILL NOT CONVERT EXACTLY ON MACHINES 01796 C WITH BASE OTHER THAN TWO, FOUR OR SIXTEEN. 01797 C THIS ROUTINE IS NOT CALLED BY ANY OTHER ROUTINE IN MP, 01798 C SO MAY BE OMITTED IF DOUBLE-PRECISION IS NOT AVAILABLE. 01799 DOUBLE PRECISION DB, DJ, DX, DBLE 01800 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01801 INTEGER B, T, R, Z(1), RS, RE, TP 01802 C CHECK LEGALITY OF B, T, M, MXR AND LUN 01803 CALL MPCHK (1, 4) 01804 I2 = T + 4 01805 C CHECK SIGN 01806 IF (DX) 20, 10, 30 01807 C IF DX = 0D0 RETURN 0 01808 10 Z(1) = 0 01809 RETURN 01810 C DX .LT. 0D0 01811 20 RS = -1 01812 DJ = -DX 01813 GO TO 40 01814 C DX .GT. 0D0 01815 30 RS = 1 01816 DJ = DX 01817 40 IE = 0 01818 50 IF (DJ.LT.1D0) GO TO 60 01819 C INCREASE IE AND DIVIDE DJ BY 16. 01820 IE = IE + 1 01821 DJ = 0.0625D0*DJ 01822 GO TO 50 01823 60 IF (DJ.GE.0.0625D0) GO TO 70 01824 IE = IE - 1 01825 DJ = 16D0*DJ 01826 GO TO 60 01827 C NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16 01828 C SET EXPONENT TO 0 01829 70 RE = 0 01830 C DB = DFLOAT(B) IS NOT ANSI STANDARD SO USE FLOAT AND DBLE 01831 DB = DBLE(FLOAT(B)) 01832 C CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) 01833 DO 80 I = 1, I2 01834 DJ = DB*DJ 01835 R(I) = IDINT(DJ) 01836 80 DJ = DJ - DBLE(FLOAT(R(I))) 01837 C NORMALIZE RESULT 01838 CALL MPNZR (RS, RE, Z, 0) 01839 IB = MAX0(7*B*B, 32767)/16 01840 TP = 1 01841 C NOW MULTIPLY BY 16**IE 01842 IF (IE) 90, 130, 110 01843 90 K = -IE 01844 DO 100 I = 1, K 01845 TP = 16*TP 01846 IF ((TP.LE.IB).AND.(TP.NE.B).AND.(I.LT.K)) GO TO 100 01847 CALL MPDIVI (Z, TP, Z) 01848 TP = 1 01849 100 CONTINUE 01850 RETURN 01851 110 DO 120 I = 1, IE 01852 TP = 16*TP 01853 IF ((TP.LE.IB).AND.(TP.NE.B).AND.(I.LT.IE)) GO TO 120 01854 CALL MPMULI (Z, TP, Z) 01855 TP = 1 01856 120 CONTINUE 01857 130 RETURN 01858 END 01859 C $$ ****** MPCHK ****** 01860 SUBROUTINE MPCHK (I, J) 01861 C MODIFIED FOR USE WITH BLAS. 01862 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01863 C CHECKS LEGALITY OF B, T, M, MXR AND LUN WHICH SHOULD BE SET 01864 C IN COMMON. 01865 C THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT 01866 C MXR .GE. (I*T + J) 01867 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01868 INTEGER B, T, R 01869 C FIRST CHECK THAT LUN IN RANGE 1 TO 99, IF NOT PRINT ERROR 01870 C MESSAGE ON LOGICAL UNIT 6. 01871 IF ((0.LT.LUN).AND.(LUN.LT.100)) GO TO 20 01872 WRITE (6, 10) LUN 01873 10 FORMAT (10H *** LUN =, I10, 26H ILLEGAL IN CALL TO MPCHK,, 01874 $ 49H PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***) 01875 LUN = 6 01876 CALL MPERR 01877 C NOW CHECK LEGALITY OF B, T AND M 01878 20 IF (B.GT.1) GO TO 40 01879 WRITE (LUN, 30) B 01880 30 FORMAT (8H *** B =, I10, 26H ILLEGAL IN CALL TO MPCHK,/ 01881 $ 49H PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***) 01882 CALL MPERR 01883 40 IF (T.GT.1) GO TO 60 01884 WRITE (LUN, 50) T 01885 50 FORMAT (8H *** T =, I10, 26H ILLEGAL IN CALL TO MPCHK,/ 01886 $ 49H PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***) 01887 CALL MPERR 01888 60 IF (M.GT.T) GO TO 80 01889 WRITE (LUN, 70) 01890 70 FORMAT (31H *** M .LE. T IN CALL TO MPCHK,/ 01891 $ 49H PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***) 01892 CALL MPERR 01893 C 8*B*B-1 SHOULD BE REPRESENTABLE, IF NOT WILL OVERFLOW 01894 C AND MAY BECOME NEGATIVE, SO CHECK FOR THIS 01895 80 IB = 4*B*B - 1 01896 IF ((IB.GT.0).AND.((2*IB+1).GT.0)) GO TO 100 01897 WRITE (LUN, 90) 01898 90 FORMAT (37H *** B TOO LARGE IN CALL TO MPCHK ***) 01899 CALL MPERR 01900 C CHECK THAT SPACE IN COMMON IS SUFFICIENT 01901 100 MX = I*T + J 01902 IF (MXR.GE.MX) RETURN 01903 C HERE COMMON IS TOO SMALL, SO GIVE ERROR MESSAGE. 01904 WRITE (LUN, 110) I, J, MX, MXR, T 01905 110 FORMAT (51H *** MXR TOO SMALL OR NOT SET TO DIM(R) BEFORE CALL, 01906 $ 21H TO AN MP ROUTINE *** / 01907 $ 27H *** MXR SHOULD BE AT LEAST, I3, 4H*T +, I4, 2H =, I6, 5H *** 01908 $ / 19H *** ACTUALLY MXR =, I10, 9H, AND T =, I10, 5H ***) 01909 CALL MPERR 01910 RETURN 01911 END 01912 C $$ ****** MPCMD ****** 01913 SUBROUTINE MPCMD (X, DZ) 01914 C MODIFIED FOR USE WITH BLAS. 01915 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01916 C CONVERTS MULTIPLE-PRECISION X TO DOUBLE-PRECISION DZ. 01917 C ASSUMES X IS IN ALLOWABLE RANGE FOR DOUBLE-PRECISION 01918 C NUMBERS. THERE IS SOME LOSS OF ACCURACY IF THE 01919 C EXPONENT IS LARGE. 01920 DOUBLE PRECISION DB, DZ, DZ2, DBLE, DLOG, DABS 01921 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01922 INTEGER B, T, R, X(1), TM 01923 C CHECK LEGALITY OF B, T, M, MXR AND LUN 01924 CALL MPCHK (1, 4) 01925 DZ = 0D0 01926 IF (X(1).EQ.0) RETURN 01927 C DB = DFLOAT(B) IS NOT ANSI STANDARD, SO USE FLOAT AND DBLE 01928 DB = DBLE(FLOAT(B)) 01929 DO 10 I = 1, T 01930 DZ = DB*DZ + DBLE(FLOAT(X(I+2))) 01931 TM = I 01932 C CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED 01933 DZ2 = DZ + 1D0 01934 C TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20, 01935 C FOR EXAMPLE ON CYBER 76. 01936 IF ((DZ2-DZ).LE.0D0) GO TO 20 01937 10 CONTINUE 01938 C NOW ALLOW FOR EXPONENT 01939 20 DZ = DZ*(DB**(X(2)-TM)) 01940 C CHECK REASONABLENESS OF RESULT. 01941 IF (DZ.LE.0D0) GO TO 30 01942 C LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG 01943 IF (DABS(DBLE(FLOAT(X(2)))-(DLOG(DZ)/ 01944 $ DLOG(DBLE(FLOAT(B)))+0.5D0)).GT.0.6D0) GO TO 30 01945 IF (X(1).LT.0) DZ = -DZ 01946 RETURN 01947 C FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL - 01948 C TRY USING MPCMDE INSTEAD. 01949 30 WRITE (LUN, 40) 01950 40 FORMAT (48H *** FLOATING-POINT OVER/UNDER-FLOW IN MPCMD ***) 01951 CALL MPERR 01952 RETURN 01953 END 01954 C $$ ****** MPDIVI ****** 01955 SUBROUTINE MPDIVI (X, IY, Z) 01956 C MODIFIED FOR USE WITH BLAS. 01957 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 01958 C DIVIDES MP X BY THE SINGLE-PRECISION INTEGER IY GIVING MP Z. 01959 C THIS IS MUCH FASTER THAN DIVISION BY AN MP NUMBER. 01960 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 01961 INTEGER B, T, R, X(1), Z(1), RS, RE, R1, C, C2, B2 01962 RS = X(1) 01963 J = IY 01964 IF (J) 30, 10, 40 01965 10 WRITE (LUN, 20) 01966 20 FORMAT (53H *** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI ***) 01967 GO TO 230 01968 30 J = -J 01969 RS = -RS 01970 40 RE = X(2) 01971 C CHECK FOR ZERO DIVIDEND 01972 IF (RS.EQ.0) GO TO 120 01973 C CHECK FOR DIVISION BY B 01974 IF (J.NE.B) GO TO 50 01975 CALL MPSTR (X, Z) 01976 IF (RE.LE.(-M)) GO TO 240 01977 Z(1) = RS 01978 Z(2) = RE - 1 01979 RETURN 01980 C CHECK FOR DIVISION BY 1 OR -1 01981 50 IF (J.NE.1) GO TO 60 01982 CALL MPSTR (X, Z) 01983 Z(1) = RS 01984 RETURN 01985 60 C = 0 01986 I2 = T + 4 01987 I = 0 01988 C IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE 01989 C LONG DIVISION. ASSUME AT LEAST 16-BIT WORD. 01990 B2 = MAX0 (8*B, 32767/B) 01991 IF (J.GE.B2) GO TO 130 01992 C LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT 01993 70 I = I + 1 01994 C = B*C 01995 IF (I.LE.T) C = C + X(I+2) 01996 R1 = C/J 01997 IF (R1) 210, 70, 80 01998 C ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT 01999 80 RE = RE + 1 - I 02000 R(1) = R1 02001 C = B*(C - J*R1) 02002 KH = 2 02003 IF (I.GE.T) GO TO 100 02004 KH = 1 + T - I 02005 DO 90 K = 2, KH 02006 I = I + 1 02007 C = C + X(I+2) 02008 R(K) = C/J 02009 90 C = B*(C - J*R(K)) 02010 IF (C.LT.0) GO TO 210 02011 KH = KH + 1 02012 100 DO 110 K = KH, I2 02013 R(K) = C/J 02014 110 C = B*(C - J*R(K)) 02015 IF (C.LT.0) GO TO 210 02016 C NORMALIZE AND ROUND RESULT 02017 120 CALL MPNZR (RS, RE, Z, 0) 02018 RETURN 02019 C HERE NEED SIMULATED DOUBLE-PRECISION DIVISION 02020 130 C2 = 0 02021 J1 = J/B 02022 J2 = J - J1*B 02023 J11 = J1 + 1 02024 C LOOK FOR FIRST NONZERO DIGIT 02025 140 I = I + 1 02026 C = B*C + C2 02027 C2 = 0 02028 IF (I.LE.T) C2 = X(I+2) 02029 IF (C-J1) 140, 150, 160 02030 150 IF (C2.LT.J2) GO TO 140 02031 C COMPUTE T+4 QUOTIENT DIGITS 02032 160 RE = RE + 1 - I 02033 K = 1 02034 GO TO 180 02035 C MAIN LOOP FOR LARGE ABS(IY) CASE 02036 170 K = K + 1 02037 IF (K.GT.I2) GO TO 120 02038 I = I + 1 02039 C GET APPROXIMATE QUOTIENT FIRST 02040 180 IR = C/J11 02041 C NOW REDUCE SO OVERFLOW DOES NOT OCCUR 02042 IQ = C - IR*J1 02043 IF (IQ.LT.B2) GO TO 190 02044 C HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR 02045 IR = IR + 1 02046 IQ = IQ - J1 02047 190 IQ = IQ*B - IR*J2 02048 IF (IQ.GE.0) GO TO 200 02049 C HERE IQ NEGATIVE SO IR WAS TOO LARGE 02050 IR = IR - 1 02051 IQ = IQ + J 02052 200 IF (I.LE.T) IQ = IQ + X(I+2) 02053 IQJ = IQ/J 02054 C R(K) = QUOTIENT, C = REMAINDER 02055 R(K) = IQJ + IR 02056 C = IQ - J*IQJ 02057 IF (C.GE.0) GO TO 170 02058 C CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED 02059 210 CALL MPCHK (1, 4) 02060 WRITE (LUN, 220) 02061 220 FORMAT (48H *** INTEGER OVERFLOW IN MPDIVI, B TOO LARGE ***) 02062 230 CALL MPERR 02063 Z(1) = 0 02064 RETURN 02065 C UNDERFLOW HERE 02066 240 CALL MPUNFL(Z) 02067 RETURN 02068 END 02069 C $$ ****** MPERR ****** 02070 SUBROUTINE MPERR 02071 C MODIFIED FOR USE WITH BLAS. 02072 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02073 C THIS ROUTINE IS CALLED WHEN A FATAL ERROR CONDITION IS 02074 C ENCOUNTERED, AND AFTER A MESSAGE HAS BEEN WRITTEN ON 02075 C LOGICAL UNIT LUN. 02076 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02077 INTEGER B, T, R 02078 WRITE (LUN, 10) 02079 10 FORMAT (42H *** EXECUTION TERMINATED BY CALL TO MPERR, 02080 $ 25H IN MP VERSION 770217 ***) 02081 C AT PRESENT JUST STOP, BUT COULD DUMP B, T, ETC. HERE. 02082 C ACTION COULD EASILY BE CONTROLLED BY A FLAG IN LABELLED COMMON. 02083 C ANSI VERSION USES STOP, UNIVAC 1108 VERSION USES 02084 C RETURN 0 IN ORDER TO GIVE A TRACE-BACK. 02085 C FOR DEBUGGING PURPOSES IT MAY BE USEFUL SIMPLY TO 02086 C RETURN HERE. MOST MP ROUTINES RETURN WITH RESULT 02087 C ZERO AFTER CALLING MPERR. 02088 STOP 02089 END 02090 C $$ ****** MPMAXR ****** 02091 SUBROUTINE MPMAXR (X) 02092 C MODIFIED FOR USE WITH BLAS. 02093 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02094 C SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER 02095 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02096 INTEGER B, T, R, X(1) 02097 C CHECK LEGALITY OF B, T, M, MXR AND LUN 02098 CALL MPCHK (1, 4) 02099 IT = B - 1 02100 C SET FRACTION DIGITS TO B-1 02101 DO 10 I = 1, T 02102 10 X(I+2) = IT 02103 C SET SIGN AND EXPONENT 02104 X(1) = 1 02105 X(2) = M 02106 RETURN 02107 END 02108 C $$ ****** MPMUL ****** 02109 SUBROUTINE MPMUL (X, Y, Z) 02110 C MODIFIED FOR USE WITH BLAS. 02111 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02112 C MULTIPLIES X AND Y, RETURNING RESULT IN Z, FOR MP X, Y AND Z. 02113 C THE SIMPLE O(T**2) ALGORITHM IS USED, WITH 02114 C FOUR GUARD DIGITS AND R*-ROUNDING. 02115 C ADVANTAGE IS TAKEN OF ZERO DIGITS IN X, BUT NOT IN Y. 02116 C ASYMPTOTICALLY FASTER ALGORITHMS ARE KNOWN (SEE KNUTH, 02117 C VOL. 2), BUT ARE DIFFICULT TO IMPLEMENT IN FORTRAN IN AN 02118 C EFFICIENT AND MACHINE-INDEPENDENT MANNER. 02119 C IN COMMENTS TO OTHER MP ROUTINES, M(T) IS THE TIME 02120 C TO PERFORM T-DIGIT MP MULTIPLICATION. THUS 02121 C M(T) = O(T**2) WITH THE PRESENT VERSION OF MPMUL, 02122 C BUT M(T) = O(T.LOG(T).LOG(LOG(T))) IS THEORETICALLY POSSIBLE. 02123 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02124 INTEGER B, T, R, X(1), Y(1), Z(1), RS, RE, XI, C, RI 02125 C CHECK LEGALITY OF B, T, M, MXR AND LUN 02126 CALL MPCHK (1, 4) 02127 I2 = T + 4 02128 I2P = I2 + 1 02129 C FORM SIGN OF PRODUCT 02130 RS = X(1)*Y(1) 02131 IF (RS.NE.0) GO TO 10 02132 C SET RESULT TO ZERO 02133 Z(1) = 0 02134 RETURN 02135 C FORM EXPONENT OF PRODUCT 02136 10 RE = X(2) + Y(2) 02137 C CLEAR ACCUMULATOR 02138 DO 20 I = 1, I2 02139 20 R(I) = 0 02140 C PERFORM MULTIPLICATION 02141 C = 8 02142 DO 40 I = 1, T 02143 XI = X(I+2) 02144 C FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST 02145 IF (XI.EQ.0) GO TO 40 02146 CALL MPMLP (R(I+1), Y(3), XI, MIN0 (T, I2 - I)) 02147 C = C - 1 02148 IF (C.GT.0) GO TO 40 02149 C CHECK FOR LEGAL BASE B DIGIT 02150 IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 90 02151 C PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME, 02152 C FASTER THAN DOING IT EVERY TIME. 02153 DO 30 J = 1, I2 02154 J1 = I2P - J 02155 RI = R(J1) + C 02156 IF (RI.LT.0) GO TO 70 02157 C = RI/B 02158 30 R(J1) = RI - B*C 02159 IF (C.NE.0) GO TO 90 02160 C = 8 02161 40 CONTINUE 02162 IF (C.EQ.8) GO TO 60 02163 IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 90 02164 C = 0 02165 DO 50 J = 1, I2 02166 J1 = I2P - J 02167 RI = R(J1) + C 02168 IF (RI.LT.0) GO TO 70 02169 C = RI/B 02170 50 R(J1) = RI - B*C 02171 IF (C.NE.0) GO TO 90 02172 C NORMALIZE AND ROUND RESULT 02173 60 CALL MPNZR (RS, RE, Z, 0) 02174 RETURN 02175 70 WRITE (LUN, 80) 02176 80 FORMAT (47H *** INTEGER OVERFLOW IN MPMUL, B TOO LARGE ***) 02177 GO TO 110 02178 90 WRITE (LUN, 100) 02179 100 FORMAT (43H *** ILLEGAL BASE B DIGIT IN CALL TO MPMUL,, 02180 $ 33H POSSIBLE OVERWRITING PROBLEM ***) 02181 110 CALL MPERR 02182 Z(1) = 0 02183 RETURN 02184 END 02185 C $$ ****** MPMUL2 ****** 02186 SUBROUTINE MPMUL2 (X, IY, Z, TRUNC) 02187 C MODIFIED FOR USE WITH BLAS. 02188 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02189 C MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z. 02190 C MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER 02191 C EVEN IF SOME DIGITS ARE GREATER THAN B-1. 02192 C RESULT IS ROUNDED IF TRUNC.EQ.0, OTHERWISE TRUNCATED. 02193 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02194 INTEGER B, T, R, X(1), Z(1), TRUNC, RE, RS 02195 INTEGER C, C1, C2, RI, T1, T3, T4 02196 RS = X(1) 02197 IF (RS.EQ.0) GO TO 10 02198 J = IY 02199 IF (J) 20, 10, 50 02200 C RESULT ZERO 02201 10 Z(1) = 0 02202 RETURN 02203 20 J = -J 02204 RS = -RS 02205 C CHECK FOR MULTIPLICATION BY B 02206 IF (J.NE.B) GO TO 50 02207 IF (X(2).LT.M) GO TO 40 02208 CALL MPCHK (1, 4) 02209 WRITE (LUN, 30) 02210 30 FORMAT (36H *** OVERFLOW OCCURRED IN MPMUL2 ***) 02211 CALL MPOVFL (Z) 02212 RETURN 02213 40 CALL MPSTR (X, Z) 02214 Z(1) = RS 02215 Z(2) = X(2) + 1 02216 RETURN 02217 C SET EXPONENT TO EXPONENT(X) + 4 02218 50 RE = X(2) + 4 02219 C FORM PRODUCT IN ACCUMULATOR 02220 C = 0 02221 T1 = T + 1 02222 T3 = T + 3 02223 T4 = T + 4 02224 C IF J*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE 02225 C DOUBLE-PRECISION MULTIPLICATION. 02226 IF (J.GE.MAX0(8*B, 32767/B)) GO TO 110 02227 DO 60 IJ = 1, T 02228 I = T1 - IJ 02229 RI = J*X(I+2) + C 02230 C = RI/B 02231 60 R(I+4) = RI - B*C 02232 C CHECK FOR INTEGER OVERFLOW 02233 IF (RI.LT.0) GO TO 130 02234 C HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY 02235 DO 70 IJ = 1, 4 02236 I = 5 - IJ 02237 RI = C 02238 C = RI/B 02239 70 R(I) = RI - B*C 02240 IF (C.EQ.0) GO TO 100 02241 C HAVE TO SHIFT RIGHT HERE AS CARRY OFF END 02242 80 DO 90 IJ = 1, T3 02243 I = T4 - IJ 02244 90 R(I+1) = R(I) 02245 RI = C 02246 C = RI/B 02247 R(1) = RI - B*C 02248 RE = RE + 1 02249 IF (C) 130, 100, 80 02250 C NORMALIZE AND ROUND OR TRUNCATE RESULT 02251 100 CALL MPNZR (RS, RE, Z, TRUNC) 02252 RETURN 02253 C HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION 02254 110 J1 = J/B 02255 J2 = J - J1*B 02256 C FORM PRODUCT 02257 DO 120 IJ = 1, T4 02258 C1 = C/B 02259 C2 = C - B*C1 02260 I = T1 - IJ 02261 IX = 0 02262 IF (I.GT.0) IX = X(I+2) 02263 RI = J2*IX + C2 02264 IS = RI/B 02265 C = J1*IX + C1 + IS 02266 120 R(I+4) = RI - B*IS 02267 IF (C) 130, 100, 80 02268 C CAN ONLY GET HERE IF INTEGER OVERFLOW OCCURRED 02269 130 CALL MPCHK (1, 4) 02270 WRITE (LUN, 140) 02271 140 FORMAT (48H *** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***) 02272 CALL MPERR 02273 GO TO 10 02274 END 02275 C $$ ****** MPNZR ****** 02276 SUBROUTINE MPNZR (RS, RE, Z, TRUNC) 02277 C MODIFIED FOR USE WITH BLAS. 02278 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02279 C ASSUMES LONG (I.E. (T+4)-DIGIT) FRACTION IN 02280 C R, SIGN = RS, EXPONENT = RE. NORMALIZES, 02281 C AND RETURNS MP RESULT IN Z. INTEGER ARGUMENTS RS AND RE 02282 C ARE NOT PRESERVED. R*-ROUNDING IS USED IF TRUNC.EQ.0 02283 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02284 INTEGER B, T, R, Z(1), RE, RS, TRUNC, B2 02285 I2 = T + 4 02286 IF (RS.NE.0) GO TO 20 02287 C STORE ZERO IN Z 02288 10 Z(1) = 0 02289 RETURN 02290 C CHECK THAT SIGN = +-1 02291 20 IF (IABS(RS).LE.1) GO TO 40 02292 WRITE (LUN, 30) 02293 30 FORMAT (43H *** SIGN NOT 0, +1 OR -1 IN CALL TO MPNZR,, 02294 $ 33H POSSIBLE OVERWRITING PROBLEM ***) 02295 CALL MPERR 02296 GO TO 10 02297 C LOOK FOR FIRST NONZERO DIGIT 02298 40 DO 50 I = 1, I2 02299 IS = I - 1 02300 IF (R(I).GT.0) GO TO 60 02301 50 CONTINUE 02302 C FRACTION ZERO 02303 GO TO 10 02304 60 IF (IS.EQ.0) GO TO 90 02305 C NORMALIZE 02306 RE = RE - IS 02307 I2M = I2 - IS 02308 DO 70 J = 1, I2M 02309 K = J + IS 02310 70 R(J) = R(K) 02311 I2P = I2M + 1 02312 DO 80 J = I2P, I2 02313 80 R(J) = 0 02314 C CHECK TO SEE IF TRUNCATION IS DESIRED 02315 90 IF (TRUNC.NE.0) GO TO 150 02316 C SEE IF ROUNDING NECESSARY 02317 C TREAT EVEN AND ODD BASES DIFFERENTLY 02318 B2 = B/2 02319 IF ((2*B2).NE.B) GO TO 130 02320 C B EVEN. ROUND IF R(T+1).GE.B2 UNLESS R(T) ODD AND ALL ZEROS 02321 C AFTER R(T+2). 02322 IF (R(T+1) - B2) 150, 100, 110 02323 100 IF (MOD(R(T),2).EQ.0) GO TO 110 02324 IF ((R(T+2)+R(T+3)+R(T+4)).EQ.0) GO TO 150 02325 C ROUND 02326 110 DO 120 J = 1, T 02327 I = T + 1 - J 02328 R(I) = R(I) + 1 02329 IF (R(I).LT.B) GO TO 150 02330 120 R(I) = 0 02331 C EXCEPTIONAL CASE, ROUNDED UP TO .10000... 02332 RE = RE + 1 02333 R(1) = 1 02334 GO TO 150 02335 C ODD BASE, ROUND IF R(T+1)... .GT. 1/2 02336 130 DO 140 I = 1, 4 02337 IT = T + I 02338 IF (R(IT) - B2) 150, 140, 110 02339 140 CONTINUE 02340 C CHECK FOR OVERFLOW 02341 150 IF (RE.LE.M) GO TO 170 02342 WRITE (LUN, 160) 02343 160 FORMAT (35H *** OVERFLOW OCCURRED IN MPNZR ***) 02344 CALL MPOVFL (Z) 02345 RETURN 02346 C CHECK FOR UNDERFLOW 02347 170 IF (RE.LT.(-M)) GO TO 190 02348 C STORE RESULT IN Z 02349 Z(1) = RS 02350 Z(2) = RE 02351 DO 180 I = 1, T 02352 180 Z(I+2) = R(I) 02353 RETURN 02354 C UNDERFLOW HERE 02355 190 CALL MPUNFL (Z) 02356 RETURN 02357 END 02358 C $$ ****** MPOVFL ****** 02359 SUBROUTINE MPOVFL (X) 02360 C MODIFIED FOR USE WITH BLAS. 02361 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02362 C CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE 02363 C EXPONENT OF MP NUMBER X WOULD EXCEED M. 02364 C AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE 02365 C AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN, 02366 C POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER 02367 C A PRESET NUMBER OF OVERFLOWS. ACTION COULD EASILY BE DETERMINED 02368 C BY A FLAG IN LABELLED COMMON. 02369 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02370 INTEGER B, T, R, X(1) 02371 C M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC. 02372 CALL MPCHK (1, 4) 02373 C SET X TO LARGEST POSSIBLE POSITIVE NUMBER 02374 CALL MPMAXR (X) 02375 WRITE (LUN, 10) 02376 10 FORMAT (45H *** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***) 02377 C TERMINATE EXECUTION BY CALLING MPERR 02378 CALL MPERR 02379 RETURN 02380 END 02381 C $$ ****** MPSTR ****** 02382 SUBROUTINE MPSTR (X, Y) 02383 C MODIFIED FOR USE WITH BLAS. 02384 C COMMON CHANGED TO NAMED COMMON, R GIVEN DIMENSION 12. 02385 C SETS Y = X FOR MP X AND Y. 02386 COMMON /MPCOM/ B, T, M, LUN, MXR, R(12) 02387 INTEGER B, T, R, X(1), Y(1), T2 02388 C SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO) 02389 J = X(1) 02390 Y(1) = J + 1 02391 IF (J.EQ.X(1)) GO TO 10 02392 C HERE X(1) AND Y(1) MUST HAVE THE SAME ADDRESS 02393 X(1) = J 02394 RETURN 02395 C HERE X(1) AND Y(1) HAVE DIFFERENT ADDRESSES 02396 10 Y(1) = J 02397 C NO NEED TO MOVE X(2), ... IF X(1) = 0 02398 IF (J.EQ.0) RETURN 02399 T2 = T + 2 02400 DO 20 I = 2, T2 02401 20 Y(I) = X(I) 02402 RETURN 02403 END
1.7.3 
