1C 2C DRIVER FOR TESTING CMLIB ROUTINES 3C BVALUE CV FC 4C LSEI 5C 6C ONE INPUT DATA CARD IS REQUIRED 7C READ(LIN,1) KPRINT,TIMES 8C 1 FORMAT(I1,E10.0) 9C 10C KPRINT = 0 NO PRINTING 11C 1 NO PRINTING FOR PASSED TESTS, SHORT MESSAGE 12C FOR FAILED TESTS 13C 2 PRINT SHORT MESSAGE FOR PASSED TESTS, FULLER 14C INFORMATION FOR FAILED TESTS 15C 3 PRINT COMPLETE QUICK-CHECK RESULTS 16C 17C ***IMPORTANT NOTE*** 18C ALL QUICK CHECKS USE ROUTINES R2MACH AND D2MACH 19C TO SET THE ERROR TOLERANCES. 20C TIMES IS A CONSTANT MULTIPLIER THAT CAN BE USED TO SCALE THE 21C VALUES OF R1MACH AND D1MACH SO THAT 22C R2MACH(I) = R1MACH(I) * TIMES FOR I=3,4,5 23C D2MACH(I) = D1MACH(I) * TIMES FOR I=3,4,5 24C THIS MAKES IT EASILY POSSIBLE TO CHANGE THE ERROR TOLERANCES 25C USED IN THE QUICK CHECKS. 26C IF TIMES .LE. 0.0 THEN TIMES IS DEFAULTED TO 1.0 27C 28C ***END NOTE*** 29C 30 COMMON/UNIT/LUN 31 COMMON/MSG/ICNT,JTEST(38) 32 COMMON/XXMULT/TIMES 33 LUN=I1MACH(2) 34 LIN=I1MACH(1) 35 ITEST=1 36C 37C READ KPRINT,TIMES PARAMETERS FROM DATA CARD.. 38C 39 READ(LIN,1) KPRINT,TIMES 401 FORMAT(I1,E10.0) 41 IF(TIMES.LE.0.) TIMES=1. 42 CALL XSETUN(LUN) 43 CALL XSETF(1) 44 CALL XERMAX(1000) 45 IF(KPRINT.LE.1) CALL XSETF(0) 46C TEST FC 47 CALL FCQX(KPRINT,IPASS) 48 ITEST=ITEST*IPASS 49C TEST LSEI 50 CALL LSEIQX(KPRINT,IPASS) 51 ITEST=ITEST*IPASS 52C 53 IF(KPRINT.GE.1.AND.ITEST.NE.1) WRITE(LUN,2) 542 FORMAT(/,' ***** WARNING -- AT LEAST ONE TEST FOR SUBLIBRARY FC FA 55 1ILED ***** ') 56 IF(KPRINT.GE.1.AND.ITEST.EQ.1) WRITE(LUN,3) 573 FORMAT(/,' ----- SUBLIBRARY FC PASSED ALL TESTS ----- ') 58 END 59 DOUBLE PRECISION FUNCTION D2MACH(I) 60 DOUBLE PRECISION D1MACH 61 COMMON/XXMULT/TIMES 62 D2MACH=D1MACH(I) 63 IF(I.EQ.1.OR. I.EQ.2) RETURN 64 D2MACH = D2MACH * DBLE(TIMES) 65 RETURN 66 END 67 SUBROUTINE FCQX (KPRINT,IPASS) 68C 69C*********************************FCQX********************************** 70C 71C QUICK CHECK SUBPROGRAM FOR THE SUBROUTINE FC. 72C DATE AUGUST, 1978. SANDIA LABS, ALBUQUERQUE. 73C LAST UPDATE -- JAN 1980 BY LEE WALTON. 74C MAR 1988 BY RON BOISVERT (FIXED DIAGNOSTIC CHECKS) 75C JUN 1992 BY RON BOISVERT (FIXED ENDPOINT OF TEST2) 76C AUTHOR R. J. HANSON 77C 78C FIT DISCRETE DATA BY AN S-SHAPED 79C CURVE. EVALUATE THE FITTED CURVE, 80C ITS FIRST TWO DERIVATIVES, AND PROBABLE 81C ERROR CURVE. 82C 83C USE SUBPROGRAM FC( ) TO OBTAIN 84C THE CONSTRAINED CUBIC B-SPLINE 85C REPRESENTATION OF THE CURVE. 86C 87C THE VALUES OF THE COEFFICIENTS OF THE B-SPLINE AS COMPUTED 88C BY FC( ) AND THE VALUES OF THE FITTED CURVE AS COMPUTED BY BVALUE 89C IN THE DE BOOR PACKAGE ARE TESTED FOR ACCURACY WITH THE EXPECTED 90C VALUES. SEE EXAMPLE PROGRAM SAND78-1291, PP.22-27. 91C 92C*********************************************************************** 93C 94C THE DIMENSIONS IN THE FOLLOWING ARRAYS ARE AS SMALL 95C AS POSSIBLE FOR THE PROBLEM BEING SOLVED. 96 DIMENSION XDATA(09),YDATA(09),SDDATA(09) 97 DIMENSION BKPT(13),XCONST(11),YCONST(11),NDERIV(11),COEFF(09) 98 DIMENSION V(51,5) 99 DIMENSION W(529),IW(30) 100 DIMENSION CHECK(51),COEFCK(9),ITEST(38) 101 COMMON /MSG/ ICNT,ITEST 102C OUTPUT UNIT TO BE USED 103 COMMON /UNIT/ NOUT 104C 105 DATA XDATA(1),XDATA(2),XDATA(3),XDATA(4),XDATA(5), 106 1 XDATA(6),XDATA(7),XDATA(8),XDATA(9) 107 2 /0.15,0.27,0.33,0.40,0.43,0.47,0.53,0.58,0.63/ 108 DATA YDATA(1),YDATA(2),YDATA(3),YDATA(4),YDATA(5), 109 1 YDATA(6),YDATA(7),YDATA(8),YDATA(9) 110 2 /0.025,0.05,0.13,0.27,0.37,0.47,0.64,0.77,0.87/ 111 DATA SDDATA(1) /0.015 /,NDATA/09/,NORD/04/,NBKPT/13/,LAST/10/ 112 DATA BKPT(1),BKPT(2),BKPT(3),BKPT(4),BKPT(5), 113 1 BKPT(6),BKPT(7),BKPT(8),BKPT(9),BKPT(10), 114 2 BKPT(11),BKPT(12),BKPT(13) 115 3 /-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8,0.9,1.0,1.1,1.2,1.3/ 116C 117C STORE THE DATA TO BE USED TO CHECK THE ACCURACY OF THE 118C COMPUTED RESULTS. SEE SAND78-1291, P.26. 119C 120 DATA COEFCK(1),COEFCK(2),COEFCK(3),COEFCK(4),COEFCK(5), 121 1 COEFCK(6),COEFCK(7),COEFCK(8),COEFCK(9)/ 1.186380846E-13, 122 2 -2.826166426E-14, -4.333929094E-15, 1.722113311E-01, 123 3 9.421965984E-01, 9.684708719E-01, 9.894902905E-01, 124 4 1.005254855E+00, 9.894902905E-01/ 125 DATA CHECK(1), CHECK(2), CHECK(3), CHECK(4), CHECK(5), 126 1 CHECK(6), CHECK(7), CHECK(8), CHECK(9)/ 127 2 2.095830752E-16, 2.870188850E-05, 2.296151081E-04, 128 3 7.749509897E-04, 1.836920865E-03, 3.587736064E-03, 129 4 6.199607918E-03, 9.844747759E-03, 1.469536692E-02/ 130 DATA CHECK(10), CHECK(11), CHECK(12), CHECK(13), CHECK(14), 131 1 CHECK(15), CHECK(16), CHECK(17), CHECK(18)/ 132 2 2.092367672E-02, 2.870188851E-02, 3.824443882E-02, 133 3 4.993466504E-02, 6.419812979E-02, 8.146039566E-02, 134 4 1.021470253E-01, 1.266835812E-01, 1.554956261E-01/ 135 DATA CHECK(19), CHECK(20), CHECK(21), CHECK(22), CHECK(23), 136 1 CHECK(24), CHECK(25), CHECK(26), CHECK(27)/ 137 2 1.890087225E-01, 2.276484331E-01, 2.718403204E-01, 138 3 3.217163150E-01, 3.762338189E-01, 4.340566020E-01, 139 4 4.938484342E-01, 5.542730855E-01, 6.139943258E-01/ 140 DATA CHECK(28), CHECK(29), CHECK(30), CHECK(31), CHECK(32), 141 1 CHECK(33), CHECK(34), CHECK(35), CHECK(36)/ 142 2 6.716759250E-01, 7.259816530E-01, 7.755752797E-01, 143 3 8.191205752E-01, 8.556270903E-01, 8.854875002E-01, 144 4 9.094402609E-01, 9.282238286E-01, 9.425766596E-01/ 145 DATA CHECK(37), CHECK(38), CHECK(39), CHECK(40), CHECK(41), 146 1 CHECK(42), CHECK(43), CHECK(44), CHECK(45)/ 147 2 9.532372098E-01, 9.609439355E-01, 9.664352927E-01, 148 3 9.704497377E-01, 9.737257265E-01, 9.768786393E-01, 149 4 9.800315521E-01, 9.831844649E-01, 9.863373777E-01/ 150 DATA CHECK(46), CHECK(47), CHECK(48), CHECK(49), CHECK(50), 151 1 CHECK(51)/ 9.894902905E-01, 9.926011645E-01, 152 2 9.954598055E-01, 9.978139804E-01, 9.994114563E-01, 153 3 1.000000000E+00/ 154C 155C BROADCAST SDDATA(1) VALUE TO ALL OF SDDATA(*). 156 CALL SCOPY(NDATA,SDDATA,0,SDDATA,1) 157 ZERO=0. 158 ONE=1. 159 NDEG=NORD-1 160C 161C WRITE THE VARIOUS CONSTRAINTS FOR 162C THE FITTED CURVE. 163 NCONST=0 164 T=BKPT(NORD) 165C 166C CONSTRAIN FUNCTION TO BE ZERO AT LEFT-MOST BREAKPOINT. 167 NCONST=NCONST+1 168 XCONST(NCONST) = T 169 YCONST(NCONST) = ZERO 170 NDERIV(NCONST) = 2+4*0 171C 172C CONSTRAIN FIRST DERIVATIVE TO BE 173C NONNEGATIVE AT LEFT-MOST BREAKPOINT. 174 NCONST=NCONST+1 175 XCONST(NCONST) = T 176 YCONST(NCONST) = ZERO 177 NDERIV(NCONST) = 1+4*1 178C 179C CONSTRAIN SECOND DERIVATIVES TO BE 180C NONNEGATIVE AT LEFT SET OF BREAKPOINTS. 181 DO 10 I = 1, 3 182 L=NDEG+I 183 T=BKPT(L) 184 NCONST=NCONST+1 185 XCONST(NCONST) = T 186 YCONST(NCONST) = ZERO 187 NDERIV(NCONST) = 1+4*2 188 10 CONTINUE 189C 190C CONSTRAIN FUNCTION VALUE AT RIGHT-MOST 191C BREAKPOINT TO BE ONE. 192 NCONST=NCONST+1 193 T=BKPT(LAST) 194 XCONST(NCONST) = T 195 YCONST(NCONST) = ONE 196 NDERIV(NCONST) = 2+4*0 197C 198C CONSTRAIN SLOPE TO AGREE AT LEFT AND 199C RIGHT-MOST BREAKPOINTS. 200 NCONST=NCONST+1 201 XCONST(NCONST) = BKPT(NORD) 202 YCONST(NCONST) = BKPT(LAST) 203 NDERIV(NCONST) = 3+4*1 204C 205C CONSTRAIN SECOND DERIVATIVES TO BE 206C NONPOSITIVE AT RIGHT SET OF BREAKPOINTS. 207 DO 20 I = 1, 4 208 NCONST=NCONST+1 209 L=LAST-4+I 210 XCONST(NCONST) = BKPT(L) 211 YCONST(NCONST) = ZERO 212 NDERIV(NCONST) = 0+4*2 213 20 CONTINUE 214C 215 IF (KPRINT .GE. 2) WRITE (NOUT,1000) 216 1000 FORMAT('1TEST OF SUBROUTINE FC') 217 ICNT=1 218 IDIGIT=-4 219C 220 IF (KPRINT .LT. 3) GO TO 30 221 CALL SVOUT(NBKPT,BKPT,'(''1ARRAY OF KNOTS.'')',IDIGIT) 222 CALL SVOUT(NDATA,XDATA,'(//'' INDEP. VAR. VALUES'')',IDIGIT) 223 CALL SVOUT(NDATA,YDATA,'('' DEPEND. VAR. VALUES'')',IDIGIT) 224 CALL SVOUT(NDATA,SDDATA,'('' DEPEND. VAR. UNCERTAINTY'')',IDIGIT) 225C 226 CALL SVOUT(NCONST,XCONST,'(''0INDEP. VAR. CONST. VALS.'')',IDIGIT) 227 CALL SVOUT(NCONST,YCONST,'('' CONST. VALUES'')',IDIGIT) 228 CALL IVOUT(NCONST,NDERIV,'('' CONST. INDICATOR'')',IDIGIT) 229C 230 30 CONTINUE 231C 232C DECLARE AMOUNT OF WORKING STORAGE ALLOCATED TO FC( ). 233 IW(1) = 529 234 IW(2) = 30 235C 236C SET MODE TO INDICATE A NEW PROBLEM 237C AND REQUEST THE VARIANCE FUNCTION. 238 MODE=2 239C 240C OBTAIN THE COEFFICIENTS OF THE B-SPLINE. 241 CALL FC(NDATA,XDATA,YDATA,SDDATA, 242 1 NORD,NBKPT,BKPT, 243 2 NCONST,XCONST,YCONST,NDERIV, 244 3 MODE, 245 4 COEFF, 246 5 W,IW) 247C 248C CHECK COEFFICIENTS 249C 250 TOL = 10.E0*SQRT(R2MACH(3)) 251 DO 40 I = 1, NDATA 252 DIFF=ABS(COEFF(I)-COEFCK(I)) 253 IF (.NOT.(DIFF .LE. TOL)) GO TO 50 254 40 CONTINUE 255 ITEST(ICNT) = 1 256 IF (KPRINT .GE. 3) WRITE (NOUT,1001) 257 1001 FORMAT('0FC( ) PASSED TEST 1') 258 GO TO 60 259C 260 50 CONTINUE 261 ITEST(ICNT) = 0 262 IF (KPRINT .GE. 2) WRITE (NOUT,1002) 263 1002 FORMAT('0FC( ) FAILED TEST 1') 264C 265 60 CONTINUE 266 K = ITEST(ICNT) 267 IF (KPRINT .EQ. 2 .AND. K .NE. 0) GO TO 70 268 IF (KPRINT .LT. 2) GO TO 70 269 CALL SVOUT(NDATA,COEFCK,'(/'' PREDICTED COEFFICIENTS OF THE B-SPLI 270 1NE FROM SAMPLE'')',IDIGIT) 271 CALL SVOUT(NDATA,COEFF,'(/'' COEFFICIENTS OF THE B-SPLINE COMPUTED 272 1 BY FC( )'')',IDIGIT) 273C 274 70 CONTINUE 275 ICNT=ICNT+1 276C 277C COMPUTE VALUE, FIRST TWO DERIVS., AND PROBABLE UNCERTAINTY. 278 N=NBKPT-NORD 279 NVAL=51 280 XVAL=ZERO 281 DX=ONE/FLOAT(NVAL-1) 282 DO 90 I = 1, NVAL 283 IF (I .EQ. NVAL) XVAL = 1.0E0 284C 285C THE FUNCTION BVALUE( ) IS IN THE DE BOOR B-SPLINE PACKAGE. 286 DO 80 J = 1, 3 287 V(I,J+1) = BVALUE(BKPT,COEFF,N,NORD,XVAL,J-1) 288 80 CONTINUE 289 V(I,1) = XVAL 290C 291C THE VARIANCE FUNCTION CV( ) IS A COMPANION 292C SUBPROGRAM TO FC( ). 293 V(I,5) = SQRT(CV(XVAL,NDATA,NCONST,NORD,NBKPT,BKPT,W)) 294 XVAL=XVAL+DX 295 90 CONTINUE 296C 297 DO 100 I = 1, NVAL 298 DIFF=ABS(V(I,2)-CHECK(I)) 299 IF (.NOT.(DIFF .LE. TOL)) GO TO 110 300 100 CONTINUE 301 ITEST(ICNT) = 1 302 IF (KPRINT .GE. 3) WRITE (NOUT,1003) 303 1003 FORMAT('0FC( ) (AND BVALUE) PASSED TEST 2') 304 GO TO 120 305C 306 110 CONTINUE 307 ITEST(ICNT) = 0 308 IF (KPRINT .GE. 2) WRITE (NOUT,1004) 309 1004 FORMAT('0FC( ) (AND BVALUE) FAILED TEST 2') 310C 311 120 CONTINUE 312 K = ITEST(ICNT) 313 IF (KPRINT .EQ. 2 .AND. K .NE. 0) GO TO 130 314 IF (KPRINT .LT. 2) GO TO 130 315C PRINT THESE VALUES. 316 CALL SMOUT(NVAL,5,NVAL,V,'(''1'',15X,''X'',10X,''FNCN'',08X, 317 1''1ST D'',07X,''2ND D'',07X,''ERROR'')', 318 2 IDIGIT) 319 WRITE (NOUT,1005) 320 1005 FORMAT('0VALUES SHOULD CORRESPOND TO THOSE IN SAND78-1291, P. 26') 321 130 CONTINUE 322C 323C CHECK ERROR PROCESSOR 324C 325 IF (KPRINT .LT. 2) GO TO 140 326 WRITE (NOUT,1006) 327 1006 FORMAT('0DIAGNOSTIC MESSAGES (6) FOR FC') 328 CALL FC(NDATA,XDATA,YDATA,SDDATA,0,NBKPT,BKPT,NCONST,XCONST, 329 1 YCONST,NDERIV,MODE,COEFF,W,IW) 330 CALL FC(NDATA,XDATA,YDATA,SDDATA,NORD,0,BKPT,NCONST,XCONST, 331 1 YCONST,NDERIV,MODE,COEFF,W,IW) 332 CALL FC(-1,XDATA,YDATA,SDDATA,NORD,NBKPT,BKPT,NCONST,XCONST, 333 1 YCONST,NDERIV,MODE,COEFF,W,IW) 334 IW(1) = 529 335 MODE = 0 336 CALL FC(NDATA,XDATA,YDATA,SDDATA,NORD,NBKPT,BKPT,NCONST,XCONST, 337 1 YCONST,NDERIV,MODE,COEFF,W,IW) 338 IW(1) = 10 339 CALL FC(NDATA,XDATA,YDATA,SDDATA,NORD,NBKPT,BKPT,NCONST,XCONST, 340 1 YCONST,NDERIV,MODE,COEFF,W,IW) 341 IW(1) = 529 342 IW(2) = 2 343 MODE = 2 344 CALL FC(NDATA,XDATA,YDATA,SDDATA,NORD,NBKPT,BKPT,NCONST,XCONST, 345 1 YCONST,NDERIV,MODE,COEFF,W,IW) 346C 347 140 CONTINUE 348 IP=1 349 DO 150 I = 1, ICNT 350 IP=IP*ITEST(I) 351 150 CONTINUE 352 IPASS=IP 353 IF (IPASS.EQ.1 .AND. KPRINT.GE.2) WRITE (NOUT,1007) 354 IF (IPASS.EQ.0 .AND. KPRINT.GE.1) WRITE (NOUT,1008) 355 1007 FORMAT(/' ***** FC PASSED ALL TESTS *****') 356 1008 FORMAT(/' ***** FC FAILED SOME TESTS *****') 357 RETURN 358 END 359 SUBROUTINE LSEIQX(KPRINT, IPASS) 360C 361C ******************************** LSEIQX ***************************** 362C 363C QUICK CHECK SUBPROGRAM FOR THE SUBROUTINE LSEI. 364C DATE FEBRUARY 16, 1979. SANDIA LABS, ALBUQUERQUE. 365C LAST UPDATE -- JANUARY, 1980 BY LEE WALTON. 366C -- April, 1996 by Ron Boisvert. Gave IP(1), IP(2) values 367C AUTHORS R. J. HANSON, KAREN HASKELL. 368C 369C THE SAMPLE PROBLEM SOLVED IS FROM A PAPER BY J. STOER, IN 370C SIAM JOURNAL OF NUM. ANAL., JUNE 1971. 371C 372C ********************************************************************* 373C 374 DIMENSION D(11,6), IP(17), WORK(105), F(6) 375 DIMENSION X(5), H(5), SOL(5), A(6,5), G(5,5) 376 DIMENSION PRGOPT(4) 377 DIMENSION ITEST(38) 378 COMMON /MSG/ ICNT, ITEST 379C OUTPUT UNIT TO BE USED. 380 COMMON /UNIT/ NOUT 381C 382C DEFINE THE DATA ARRAYS FOR THE EXAMPLE. THE ARRAY A( ) 383C CONTAINS THE LEAST SQUARES EQUATIONS. (THERE ARE NO EQUALITY 384C CONSTRAINTS IN THIS EXAMPLE). 385 DATA A(1,1), A(1,2), A(1,3), A(1,4), A(1,5) /-74.,80.,18.,-11., 386 * -4./ 387 DATA A(2,1), A(2,2), A(2,3), A(2,4), A(2,5) /14.,-69.,21.,28.,0./ 388 DATA A(3,1), A(3,2), A(3,3), A(3,4), A(3,5) /66.,-72.,-5.,7.,1./ 389 DATA A(4,1), A(4,2), A(4,3), A(4,4), A(4,5) /-12.,66.,-30.,-23., 390 * 3./ 391 DATA A(5,1), A(5,2), A(5,3), A(5,4), A(5,5) /3.,8.,-7.,-4.,1./ 392 DATA A(6,1), A(6,2), A(6,3), A(6,4), A(6,5) /4.,-12.,4.,4.,0./ 393C 394C THE ARRAY G( ) CONTAINS THE INEQUALITY CONSTRAINT EQUATIONS, 395C WRITTEN IN THE SENSE 396C (ROW VECTOR)*(SOLUTION VECTOR) .GE. (GIVEN VALUE). 397 DATA G(1,1), G(1,2), G(1,3), G(1,4), G(1,5) /-1.,-1.,-1.,-1.,-1./ 398 DATA G(2,1), G(2,2), G(2,3), G(2,4), G(2,5) /10.,10.,-3.,5.,4./ 399 DATA G(3,1), G(3,2), G(3,3), G(3,4), G(3,5) /-8.,1.,-2.,-5.,3./ 400 DATA G(4,1), G(4,2), G(4,3), G(4,4), G(4,5) /8.,-1.,2.,5.,-3./ 401 DATA G(5,1), G(5,2), G(5,3), G(5,4), G(5,5) /-4.,-2.,3.,-5.,1./ 402C 403C DEFINE THE LEAST SQUARES RIGHT-SIDE VECTOR. 404 DATA F(1), F(2), F(3), F(4) /-5.,-9.,708.,4165./ 405 DATA F(5), F(6) /-13266.,8409./ 406C 407C DEFINE THE INEQUALITY CONSTRAINT RIGHT-SIDE VECTOR. 408 DATA H(1), H(2), H(3), H(4), H(5) /-5.,20.,-40.,11.,-30./ 409C 410C DEFINE THE VECTOR THAT IS THE KNOWN SOLUTION. 411 DATA SOL(1), SOL(2), SOL(3), SOL(4), SOL(5) /1.,2.,-1.,3.,-4./ 412C 413C DEFINE THE MATRIX DIMENSIONS, NUMBER OF LEAST SQUARES EQUATIONS, 414C NUMBER OF EQUALITY CONSTRAINTS, TOTAL NUMBER OF 415C EQUATIONS, AND NUMBER OF VARIABLES. SET ME=0 TO INDICATE 416C THERE ARE NO EQUALITY CONSTRAINTS. 417 MDD = 11 418 MDA = 6 419 MDG = 5 420 MA = 6 421 MG = 5 422 M = MA + MG 423 N = 5 424 ME = 0 425C 426 ip(1) = 105 427 ip(2) = 17 428C 429 NP1 = N + 1 430 MEP1 = ME + 1 431 MEAP1 = ME + MA + 1 432C 433C COPY THE PROBLEM MATRICES 434 DO 10 I = 1, N 435C 436C COPY THE I-TH COL OF THE INEQUALITY CONSTRAINT MATRIX INTO 437C THE WORK ARRAY. 438 CALL SCOPY(MG, G(1,I), 1, D(MEAP1,I), 1) 439C 440C COPY THE I-TH COL OF THE LEAST SQUARES MATRIX INTO THE WORK 441C ARRAY. 442 CALL SCOPY(MA, A(1,I), 1, D(MEP1,I), 1) 443C 444 10 CONTINUE 445C 446C COPY THE RIGHT-SIDE VECTORS INTO THE WORK ARRAY IN COMPATIBLE 447C ORDER. 448 CALL SCOPY(MG, H, 1, D(MEAP1,NP1), 1) 449 CALL SCOPY(MA, F, 1, D(MEP1,NP1), 1) 450C 451 ICNT = 1 452 IF (KPRINT.GE.2) WRITE (NOUT,99999) 45399999 FORMAT ('1TEST OF SUBROUTINE LSEI') 454C 455C USE DEFAULT PROGRAM OPTIONS IN LSEI( ), AND SET MATRIX- 456C VECTOR PRINTING ACCURACY PARAMETERS. 457 PRGOPT(1) = 1 458 IDIGIT = -4 459 JDIGIT = -11 460C 461C PRINT THE DATA MATRIX TO BE PASSED TO LSEI( ). 462C CALL SMOUT (MA,N+1,MDD,D,'(/'' LEAST SQUARES EQUATIONS (ROWS 1-6) 463C 1DATA VALUES (COL 6)'')',IDIGIT) 464C CALL SSMOUT (MA+1,1,M,N+1,MDD,D,'('' INEQUALITY CONSTRAINT EQUATIONS 465C 1(ROWS 7-11) CONSTRAINING VALUES (COL 6)'')',IDIGIT) 466C 467C COMPUTE RESIDUAL NORM OF KNOWN LEAST SQUARES SOLN. (TO BE 468C USED TO CHECK COMPUTED RESIDUAL NORM = RNORML.) 469 DO 20 I = 1, MA 470 WORK(I) = SDOT(N,D(I,1),MDD,SOL,1) - F(I) 471 20 CONTINUE 472 RESNRM = SNRM2(MA,WORK,1) 473C 474C CALL LSEI( ) TO GET SOLN IN X(*), LEAST SQUARES RESIDUAL IN 475C RNORML. 476 CALL LSEI(D, MDD, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, MODE, 477 * WORK, IP) 478C 479C 480C COMPUTE REL. ERROR IN PROBLEM VARIABLE SOLN. AND RESIDUAL 481C NORM COMPUTATION. 482 TNORM = SNRM2(N,SOL,1) 483 CALL SAXPY(N, -1., X, 1, SOL, 1) 484 CNORM = SNRM2(N,SOL,1) 485 RELERR = CNORM/TNORM 486 RELNRM = (RESNRM-RNORML)/RESNRM 487C 488 IF (.NOT.(RELERR.LE.100.*SQRT(R2MACH(3)))) GO TO 30 489 IF (.NOT.(RELNRM.LE.10.*R2MACH(3))) GO TO 30 490 IF (KPRINT.EQ.3) WRITE (NOUT,99998) 49199998 FORMAT ('0LSEI PASSED TEST') 492 ITEST(ICNT) = 1 493 GO TO 40 494C 495 30 CONTINUE 496 ITEST(ICNT) = 0 497 IF (KPRINT.GE.2) WRITE (NOUT,99997) 49899997 FORMAT ('0LSEI FAILED TEST') 499C 500 40 CONTINUE 501 IF (KPRINT.LT.3) GO TO 50 502C PRINT OUT KNOWN SOLUTION AND COMPUTED SOLUTION 503 CALL SVOUT(N, SOL, '('' KNOWN LEAST SQUARES SOLN'')', IDIGIT) 504 CALL SVOUT(N, X, '(/'' SOLN COMPUTED BY LSEI( ).'')', JDIGIT) 505C 506 50 CONTINUE 507 K = ITEST(ICNT) 508 IF (KPRINT.EQ.2 .AND. K.NE.0) GO TO 60 509 IF (KPRINT.LT.2) GO TO 60 510C PRINT OUT THE KNOWN AND COMPUTED RESIDUAL NORMS 511 CALL SVOUT(1, RESNRM, 512 * '(/'' RESIDUAL NORM OF KNOWN LEAST SQUARES SOLN'')', JDIGIT) 513 CALL SVOUT(1, RNORML, '(/'' RES NORM COMPUTED BY LSEI( ) '')', 514 * JDIGIT) 515C PRINT OUT THE COMPUTED SOLUTION RELATIVE ERROR 516 CALL SVOUT(1, RELERR, '(/'' COMPUTED SOLN REL. ERROR'')', IDIGIT) 517C PRINT OUT THE COMPUTED RELATIVE ERROR IN RESIDUAL NORM 518 CALL SVOUT(1, RELNRM, 519 * '(/'' COMPUTED REL. ERROR IN RESIDUAL NORM'')', IDIGIT) 520C 521 60 CONTINUE 522C 523C CHECK CALLS TO ERROR PROCESSOR 524C 525 IF (KPRINT.LT.2) GO TO 70 526 WRITE (NOUT,99996) 52799996 FORMAT (/' DIAGNOSTIC MESSAGES (3) FOR LSEI ') 528 CALL LSEI(D, 0, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, MODE, 529 * WORK, IP) 530C PRGOPT(1) = 4 531C PRGOPT(2) = 1 532C PRGOPT(3) = 1 533C PRGOPT(4) = 1 534C CALL LSEI (D,MDD,ME,MA,MG,12,PRGOPT,X,RNORME,RNORML,MODE,WORK,IP) 535 PRGOPT(1) = -1 536 CALL LSEI(D, MDD, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, MODE, 537 * WORK, IP) 538C 539 70 CONTINUE 540 IPASS = ITEST(ICNT) 541 IF (IPASS.EQ.1 .AND. KPRINT.GE.2) WRITE (NOUT,99995) 542 IF (IPASS.EQ.0 .AND. KPRINT.GE.1) WRITE (NOUT,99994) 54399995 FORMAT (/' ***** LSEI PASSED ALL TESTS *****') 54499994 FORMAT (/' ***** LSEI FAILED SOME TESTS ***** ') 545C 546 RETURN 547 END 548 REAL FUNCTION R2MACH(I) 549 COMMON/XXMULT/TIMES 550 R2MACH=R1MACH(I) 551 IF(I.EQ.1.OR. I.EQ.2) RETURN 552 R2MACH = R2MACH * TIMES 553 RETURN 554 END 555 SUBROUTINE SMOUT(M,N,LDA,A,IFMT,IDIGIT) 556C ... DECLARATION REPLACED BY RON BOISVERT ON 6-NOV-91 557C DIMENSION A(LDA,N),IFMT(1),ICOL(3) 558 DIMENSION A(LDA,N),ICOL(3) 559 CHARACTER*(*) IFMT 560C ... 561C 562C SINGLE PRECISION MATRIX OUTPUT ROUTINE. 563C 564C INPUT.. 565C 566C M,N,LDA,A(*,*) PRINT THE SINGLE PRECISION ARRAY A(I,J),I = 1,...,M, 567C J=1,...,N, ON OUTPUT UNIT LOUT=6. LDA IS THE DECLARED 568C FIRST DIMENSION OF A(*,*) AS SPECIFIED IN THE CALLING 569C PROGRAM. THE HEADING IN THE FORTRAN FORMAT STATEMENT 570C IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST STEP. 571C THE COMPONENTS A(I,J) ARE INDEXED, ON OUTPUT, IN A 572C PLEASANT FORMAT. 573C IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON 574C OUTPUT UNIT LOUT=6 WITH THE VARIABLE FORMAT FORTRAN 575C STATEMENT 576C WRITE(LOUT,IFMT). 577C IDIGIT PRINT AT LEAST IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. 578C THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,10, OR 14 579C WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF 580C PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE 581C UTILIZED TO WRITE EACH LINE OF OUTPUT OF THE ARRAY 582C A(*,*). (THIS CAN BE USED ON MOST TIME-SHARING 583C TERMINALS). IF IDIGIT.GE.0, 133 PRINTING COLUMNS ARE 584C UTILIZED. (THIS CAN BE USED ON MOST LINE PRINTERS). 585C 586C EXAMPLE.. 587C 588C PRINT AN ARRAY CALLED (SIMPLEX TABLEAU ) OF SIZE 10 BY 20 SHOWING 589C 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING 590C SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. 591C 592C DIMENSION TABLEU(20,20) 593C M = 10 594C N = 20 595C LDTABL = 20 596C IDIGIT = -6 597C CALL SMOUT(M,N,LDTABL,TABLEU,21H(16H1SIMPLEX TABLEAU),IDIGIT) 598C 599C 600C 601C AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. 602C RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. 603C DATE JULY 27,1978. 604C 605C 606 DATA ICOL(1),ICOL(2),ICOL(3)/1HC,1HO,1HL/ 607 LOUT=I1MACH(2) 608 WRITE(LOUT,IFMT) 609 IF(M.LE.0.OR.N.LE.0.OR.LDA.LE.0) RETURN 610 NDIGIT = IDIGIT 611 IF(IDIGIT.EQ.0) NDIGIT = 4 612 IF(IDIGIT.GE.0) GO TO 80 613C 614 NDIGIT = -IDIGIT 615 IF(NDIGIT.GT.4) GO TO 20 616C 617 DO 10 K1=1,N,5 618 K2 = MIN0(N,K1+4) 619 WRITE(LOUT,1000) (ICOL,I,I = K1, K2) 620 DO 10 I = 1, M 621 WRITE(LOUT,1004) I,(A(I,J),J = K1, K2) 622 10 CONTINUE 623 RETURN 624C 625 20 CONTINUE 626 IF(NDIGIT.GT.6) GO TO 40 627C 628 DO 30 K1=1,N,4 629 K2 = MIN0(N,K1+3) 630 WRITE(LOUT,1001) (ICOL,I,I = K1, K2) 631 DO 30 I = 1, M 632 WRITE(LOUT,1005) I,(A(I,J),J = K1, K2) 633 30 CONTINUE 634 RETURN 635C 636 40 CONTINUE 637 IF(NDIGIT.GT.10) GO TO 60 638C 639 DO 50 K1=1,N,3 640 K2=MIN0(N,K1+2) 641 WRITE(LOUT,1002) (ICOL,I,I = K1, K2) 642 DO 50 I = 1, M 643 WRITE(LOUT,1006) I,(A(I,J),J = K1, K2) 644 50 CONTINUE 645 RETURN 646C 647 60 CONTINUE 648 DO 70 K1=1,N,2 649 K2 = MIN0(N,K1+1) 650 WRITE(LOUT,1003) (ICOL,I,I = K1, K2) 651 DO 70 I = 1, M 652 WRITE(LOUT,1007) I,(A(I,J),J = K1, K2) 653 70 CONTINUE 654 RETURN 655C 656 80 CONTINUE 657 IF(NDIGIT.GT.4) GO TO 100 658C 659 DO 90 K1=1,N,10 660 K2 = MIN0(N,K1+9) 661 WRITE(LOUT,1000) (ICOL,I,I = K1, K2) 662 DO 90 I = 1, M 663 WRITE(LOUT,1004) I,(A(I,J),J = K1, K2) 664 90 CONTINUE 665 RETURN 666C 667 100 CONTINUE 668 IF(NDIGIT.GT.6) GO TO 120 669C 670 DO 110 K1=1,N,8 671 K2 = MIN0(N,K1+7) 672 WRITE(LOUT,1001) (ICOL,I,I = K1, K2) 673 DO 110 I = 1, M 674 WRITE(LOUT,1005) I,(A(I,J),J = K1, K2) 675 110 CONTINUE 676 RETURN 677C 678 120 CONTINUE 679 IF(NDIGIT.GT.10) GO TO 140 680C 681 DO 130 K1=1,N,6 682 K2 = MIN0(N,K1+5) 683 WRITE(LOUT,1002) (ICOL,I,I = K1, K2) 684 DO 130 I = 1, M 685 WRITE(LOUT,1006) I,(A(I,J),J = K1, K2) 686 130 CONTINUE 687 RETURN 688C 689 140 CONTINUE 690 DO 150 K1=1,N,5 691 K2 = MIN0(N,K1+4) 692 WRITE(LOUT,1003) (ICOL,I,I = K1, K2) 693 DO 150 I = 1, M 694 WRITE(LOUT,1007) I,(A(I,J),J = K1, K2) 695 150 CONTINUE 696 RETURN 697 1000 FORMAT(10X,10(4X,3A1,I4,1X)) 698 1001 FORMAT(10X,8(5X,3A1,I4,2X)) 699 1002 FORMAT(10X,6(7X,3A1,I4,4X)) 700 1003 FORMAT(10X,5(9X,3A1,I4,6X)) 701 1004 FORMAT(1X,'ROW',I4,2X,1P10E12.3) 702 1005 FORMAT(1X,'ROW',I4,2X,1P8E14.5) 703 1006 FORMAT(1X,'ROW',I4,2X,1P6E18.9) 704 1007 FORMAT(1X,'ROW',I4,2X,1P5E22.13) 705 END 706