C ALGORITHM 645 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 3, C SEPT., 1986, P. 274. C MAINLINE PROGRAM TO TEST GENERALIZED INVERSE SUBROUTINES C C J.C.NASH 1979, J.C.NASH AND C.E.GRATTON 1982 C J.C.NASH 1983, 1984 C J.C.NASH AND R.L.C.WANG 1985, 1986 C C INTEGER I, ISEED, J, K, M, MA, MOPT, N, NA, NB, NIN, NOUT LOGICAL FAIL REAL A(30,30), C(30,30), G(30,30), TA(4), TM(4), X(30,30) REAL B(30,30) REAL ALPHA, TOL INTEGER MSG(80) C C DATA YES/1HY/ C C THE ABOVE STATEMENT IS USED IN FORTRAN LEVEL 66 AND WILL C PASS THE PFORT VERIFIER. HOWEVER, SOME COMPILERS, FOR C EXAMPLE THE MICROSOFT FORTRAN COMPILER (VERSION 3.2), C REQUIRE THAT CHARACTER DATA BE DECLARED AS SUCH. THE C FOLLOWING CODE WAS USED FOR THE MICROSOFT COMPILER -- C C CHARACTER YES,TEST DATA YES/'Y'/ C C NOTE THAT TEST AND YES MUST BE DELETED FROM THE LAST C REAL VARIABLE DECLARATION ABOVE C C C I/O CHANNELS C C INCLUDE THE NEXT TWO STATEMENTS FOR MICROSOFT FORTRAN OPEN (5,FILE=' ') OPEN (6,FILE=' ') C C NIN = 5 NOUT = 6 C C READ A CHARACTER TO DECIDE IF CALLING SEQUENCE TESTS ARE TO C BE CARRIED OUT. NOTE THAT SOME FORTRAN PROCESSORS, IN C PARTICULAR WATFIV, WILL GENERATE AN EXECUTION ERROR WHICH C HALTS THE PROGRAM, SO THAT THESE TESTS SHOULD NOT BE RUN. C C THE CHARACTER MUST BE THE LETTER Y FOR THE TESTS TO BE C CARRIED OUT. ANY OTHER CHARACTERS WILL CAUSE THE TESTS TO BE C BYPASSED. (CHARACTERS APART FROM THE FIRST ARE IGNORED.) C READ (NIN, 9998) TEST IF (TEST .NE. YES) GO TO 5 C C C THE FOLLOWING SEQUENCE OF CALLS IS USED TO TEST ERROR TRAPS IN C GINVS, PTST AND ZIELKE C THIS SECTION CAN BE OMITTED WITHOUT AFFECTING COMPUTATIONS C WRITE (NOUT,9991) C C TEST TRAP FOR MATRIX ROW DIMENSION LESS THAN 1 C C TESTING (MA.LT.1) C MA = 0 NA = 10 M = 10 N = 10 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) IF (.NOT.FAIL) STOP C C TEST TRAP FOR MATRIX COLUMN DIMENSION LESS THAN 1 C C TESTING (NA.LT.1) C MA = 11 NA = 0 M = 11 N = 11 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) IF (.NOT.FAIL) STOP C C TEST TRAP FOR MATRIX ROW ORDER LESS THAN 1 C C TESTING (M.LT.1) C MA = 12 NA = 12 M = 0 N = 12 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) IF (.NOT.FAIL) STOP C C TEST TRAP FOR MATRIX COLUMN ORDER LESS THAN 1 C C TESTING (N.LT.1) C MA = 13 NA = 13 M = 13 N = 0 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) IF (.NOT.FAIL) STOP C C TEST TRAP FOR NUMBER OF ROWS GREATER THAN DIMENSION C C TESTING (M.GT.MA) C MA = 5 NA = 5 M = 6 N = 5 CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP C C TEST TRAP FOR NUMBER OF COLUMNS GREATER THAN DIMENSION C C TESTING (N.GT.NA) C MA = 5 NA = 5 M = 5 N = 7 CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP C C TESTING MA AND NA FORBIDDEN BY ZIELKE C MA = 3 NA = 3 M = 3 N = 3 MOPT = 2 CALL ZIELKE(M, N, A, MA, NA, MOPT, ALPHA, NOUT, FAIL) C C OUGHT TO FAIL, STOP IF IT HAS NOT FAILED C IF (.NOT.FAIL) STOP C C END OF SEQUENCE OF CALLS TO TEST ERROR TRAPS WRITE (NOUT,9992) C C DECIDE IF FIXED TEST OF ROUTINE PTST IS TO BE RUN. C C READ A CHARACTER TO DECIDE IF FIXED TESTS OF PTST ARE TO TO C BE CARRIED OUT. C C THE CHARACTER MUST BE THE LETTER Y FOR THE TESTS TO BE C CARRIED OUT. ANY OTHER INPUT WILL CAUSE THE TESTS TO BE C BYPASSED. (CHARACTERS APART FROM THE FIRST ARE IGNORED.) C 5 READ (NIN, 9998) TEST IF(TEST .NE. YES) GO TO 10 C C C C THE FOLLOWING TESTS ARE DESIGNED TO ENSURE THAT THE ROUTINE C PTST IS PERFORMING CORRECTLY. C FIRST SET UP THE MATRIX AND A SUPPOSED INVERSE. C C MATRIX C WRITE (NOUT, 9990) C M = 3 N = 2 A(1,1)=1. A(1,2)=2. A(2,1)=2. A(2,2)=-3. A(3,1)=0.5 A(3,2)=0.0 C C SUPPOSED INVERSE C X(1,1)=0.1 X(1,2)=0.2 X(1,3)=-0.2 X(2,1)=0.3 X(2,2)=-0.3 X(2,3)=0.9 MA = 30 NA = 30 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C C APPROXIMATE RESULTS C C TEST PTST WITH 3 BY 2 MATRICES C TO VERIFY MATRIX MULTIPLICATIONS C C A X C C 1. 2. 0.1 0.2 -0.2 C 2. -3. 0.3 -0.3 0.9 C 0.5 0. C C TEST PENROSE CONDITIONS C C A IS THE INPUT MATRIX, X IS THE INVERSE OF A C INPUT MATRIX NORM= 4.27200100E 00 C INVERSE MATRIX NORM= 1.03923000E 00 C AVERAGE DEVIATION MAXIMUM DEVIATION C AXA=TEST=A ... ACTUAL 0.891666 2.3 C NORMALIZED 0.2087232 0.5383886 C XAX=TEST=X ... ACTUAL 0.1875 0.42 C NORMALIZED 0.1804217 0.4041449 C (AX)T=TEST=AX ...ACTUAL 0.1683332 0.32 C NORMALIZED 0.3791636 0.7207864 C (XA)T=TEST=XA ...ACTUAL 0.55 0.55 C NORMALIZED 0.1238851 0.1238851 C C REVERSE CALL TO TEST ROWSIZE .LT. COLUMNSIZE CASE C WRITE (NOUT, 9989) CALL PTST(N, M, X, NA, A, MA, C, MA, TA, TM, NOUT, FAIL) C C C TEST ROW AND COLUMN VECTOR PROBLEMS C C MATRIX (COLUMN VECTOR) A C WRITE (NOUT, 9988) C A(1,1) = 3. A(2,1) = 4. C C MATRIX (ROW VECTOR) X C X(1,1) = 3. / 25. X(1,2) = 4. / 25. C M = 2 N = 1 CALL PTST(M, N, A, MA, X, NA, C, MA, TA, TM, NOUT, FAIL) C WRITE (NOUT, 9989) C C REVERSE CALL TO TEST ROWSIZE .LT. COLUMNSIZE CASE C CALL PTST(N, M, X, NA, A, MA, C, MA, TA, TM, NOUT, FAIL) C WRITE (NOUT, 9987) C C TEST TRIVIAL CASE (1 BY 1) C A(1,1) = 2. X(1,1) = 0.5 M = 1 N = 1 C CALL PTST(N, M, X, NA, A, MA, C, MA, TA, TM, NOUT, FAIL) C TEST MAIN DRIVER CODE FOR GENERALIZED INVERSE TESTING C C MA IS MAXIMUM FIRST DIMENSION OF WORKING ARRAYS C NA IS MAXIMUM SECOND DIMENSION OF WORKING ARRAYS C NB IS MAX(MA,NA) AND THE FIRST AND SECOND DIMENSION OF THE C WORKING ARRAY B. C C SIZE OF THE INPUT MATRIX SHOULD NOT EXCEED 30 BY 30 C C TOP OF CYCLE C 10 MA = 30 NA = 30 NB = 30 C C READ TWO RECORDS WHICH CONTROL WHICH TEST IS PERFORMED C C FIRST RECORD IS A TEXT MESSAGE DESCRIBING THE TEST C MESSAGE MAY BE UP TO 80 CHARACTERS IN LENGTH C C SECOND RECORD CONTAINS CONTROL INFORMATION C M AND N GIVE THE SIZE OF THE MATRIX C TO BE GENERATED BY GMATX C K GIVES ITS RANK C MOPT = 0 IF GMATX TO BE CALLED C MOPT = 1, 2, 3, 4 FOR ZIELKE MATRICES C =-1,-2,-3,-4 FOR THEIR MOORE-PENROSE INVERSES C ISEED = AN INTEGER SEED FOR USE BY THE RANDOM C NUMBER GENERATOR CALLED BY GMATX C ALPHA = A PARAMETER USED TO ADJUST THE C SINGULAR VALUES OF THE MATRIX GENERATED C BY GMATX C C FOR DETAILS OF THE CONTROLS, SEE THE COMMENTS IN THE C SUBROUTINES GMATX AND ZIELKE C C FOR THE FORMATS OF THE INPUT FIELDS, SEE FORMAT 9998 C BELOW. C C READ (NIN,9998) (MSG(I),I=1,80) WRITE (NOUT,9996) WRITE (NOUT,9997) (MSG(I),I=1,80) READ (NIN,9995) M, N, K, MOPT, ISEED, ALPHA C C INITIALIZE ERROR FLAG TO IMPLY CORRECT EXECUTION C FAIL = .FALSE. C C CHECK FOR END OF DATA C IF (.NOT.(M.EQ.0 .AND. N.EQ.0 .AND. K.EQ.0 .AND. * MOPT.EQ.0 .AND. ISEED.EQ.0 .AND. ALPHA.EQ.(0.0))) GO TO * 20 WRITE (NOUT,9993) STOP 20 WRITE (NOUT,9994) M, N, K, MOPT, ISEED, ALPHA C OTHER TESTS FOR VALID INPUTS ARE MADE IN THE SUBROUTINES C C CHECK FOR DIMENSIONS EXCEEDED C IF (M.LT.MA .AND. N.LT.NA) GO TO 30 WRITE (NOUT,9999) GO TO 10 30 IF (MOPT.EQ.0) GO TO 40 C C ZIELKE ROUTINE CALL C C NOTE THAT M AND N ARE REPLACED BY APPROPRIATE VALUES, C MOPT IS THE ZIELKE MATRIX SELECTED C ALPHA IS NEEDED IN ZIELKE AS A PARAMETER IN C FORMULAS FOR THE MATRIX ELEMENTS GENERATED. C CALL ZIELKE(M, N, A, MA, NA, MOPT, ALPHA, NOUT, FAIL) C C QUESTION - WAS AN ERROR DETECTED C IF (FAIL) GO TO 10 GO TO 50 C C USING GENERATED MATRIX C C WE SUGGEST VALUES OF ALPHA TO LIE WITHIN 0.10 AND 10.0, BUT C OTHER VALUES WILL BE TOLERATED. NOTE THAT ALPHA CANNOT BE C EQUAL TO ZERO, AND IF SO WILL BE SET TO ONE IN GMATX C 40 CALL GMATX(M, N, A, MA, NA, B, NB, ALPHA, K, NOUT, ISEED, FAIL) C C QUESTION - WAS AN ERROR DETECTED C IF (FAIL) GO TO 10 C C SUBROUTINE CALL FOR GENERALIZED INVERSE CALCULATOR C THIS CALL MUST BE CHANGED FOR EACH ROUTINE TESTED C C GINVSE IS NOT DESIGNED TO HANDLE MORE COLUMNS THAN ROWS C THAT IS, N .GT. M, SO AVOID CALL IN SUCH CASES C 50 IF (N.GT.M) GO TO 10 C C GINVSE IS ALSO UNABLE TO HANDLE COLUMN VECTORS, THAT C IS AN M BY 1 MATRIX A. C IF (N .LT. 2) GO TO 10 C C C NOTE ADDITIOMAL STATEMENTS TO GET CONFORMITY WITH TEST PROGRAM C COPY MATRIX TO AVOID OVERWRITING IT C DO 70 I=1,M DO 60 J=1,N G(I,J) = A(I,J) 60 CONTINUE 70 CONTINUE C C IBM SINGLE PRECISION C TOL = 16.0**(-5) CALL GINVSE(G, MA, NA, X, MA, M, N, TOL, C, MA, FAIL, * NOUT) C C QUESTION - WAS AN ERROR DETECTED C IF (FAIL) GO TO 10 C C END OF SPECIAL CODE FOR GINVSE, THE PROGRAM UNDER TEST C C TEST PENROSE CONDITIONS C CALL PTST(M, N, A, MA, X, MA, C, MA, TA, TM, NOUT, FAIL) IF (FAIL) GO TO 10 C C RETURN TO TOP OF CYCLE C GO TO 10 C C * * * * * * * * * * * * * * * C 9999 FORMAT (45H0**ERROR** DIMENSIONS EXCEEDED, INPUT IGNORED) 9998 FORMAT (80A1) 9997 FORMAT (56A1/24A1) 9996 FORMAT (20H0 DATA VERIFICATION) 9995 FORMAT (4I5, I10, F10.5) 9994 FORMAT (41H0 M N K MOPT ISEED ALPHA / * 1X, 4I5, I10, F10.5, 2(/)) 9993 FORMAT (5(/), 40H END OF DATA HAS BEEN REACHED (ALL INPUT, * 34H VALUES ARE CONSIDERED TO BE ZERO), 5(/)) 9992 FORMAT (38H0 END OF CALLING SEQUENCE ERROR TESTS) 9991 FORMAT (41H0 TEST TRAPS FOR CALLING SEQUENCE ERRORS / * 44H PROVOKE 11 ERROR REPORTS FOR DIMENSIONS AND, * 18H INVALID ARGUMENTS ) 9990 FORMAT (31H0TEST PTST WITH 3 BY 2 MATRICES / * 33H TO VERIFY MATRIX MULTIPLICATIONS // * 10X, 2H A, 27X, 2H X // * 8X, 9H 1. 2., 15X, 16H 0.1 0.2 -0.2 / * 8X, 9H 2. -3., 15X, 16H 0.3 -0.3 0.9 / * 8X, 9H 0.5 0. ) 9989 FORMAT (33H0SAME TEST WITH A AND X EXCHANGED ) 9988 FORMAT (38H0TEST PTST WITH ROW AND COLUMN VECTORS / * 10X, 2H A, 23X, 2H X / * 8X, 5H 3.0, 16X, 15H 0.12 0.16 / * 8X, 5H 4.0 ) 9987 FORMAT (39H0TEST WITH TRIVIAL 1 BY 1 MATRIX A = 2. ) C C * * * END OF MAIN PROGRAM * * * C END DOUBLE PRECISION FUNCTION DRAND(IX) C C PORTABLE FORTRAN RANDOM NUMBER GENERATOR C USING THE RECURSION C IX = IX * A MOD P C C FROM LINUS SHRAGE, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL 5, NO 2, P 134 FF, JUNE 1979 C C CALLING SEQUENCE -- C C DPPRN = DRAND ( IX ) C C IX = SEED FOR NEXT MEMBER OF PSEUDO-RANDOM SEQUENCE C SHOULD NOT BE ALTERED BETWEEN CALLS TO DRAND C NOTE -- IX IS A DOUBLE PRECISION VARIABLE C IX IS ALTERED ON EACH CALL TO DRAND C DPPRN = DOUBLE PRECISION PSEUDO RANDOM NUMBER IN THE C UNIT INTERVAL WHICH IS PRODUCED BY THIS C FUNCTION C C UNIFORM PSEUDO RANDOM NUMBER ON UNIT INTERVAL IS C RETURNED VIA DRAND (DOUBLE PRECISION) C DOUBLE PRECISION A, P, IX, B15, B16, XHI, XALO, LEFTLO, * FHI, K C C 7**15, 2**15, 2**16, 2**31 - 1 C DATA A /16807.D0/, B15 /32768.D0/, B16 /65536.D0/, P * /2147483647.D0/ C C GET 15 HI ORDER BITS OF IX C XHI = IX/B16 XHI = XHI - DMOD(XHI,1.0D0) C C GET 16 LO BITS OF IX AND FORM LO PRODUCT C XALO = (IX-XHI*B16)*A C C GET 15 HI ORDER BITS OF LO PRODUCT C LEFTLO = XALO/B16 LEFTLO = LEFTLO - DMOD(LEFTLO,1.0D0) C C FORM THE 31 HIGHEST BITS OF FULL PRODUCT C FHI = XHI*A + LEFTLO C C GET OVERFLOW PAST 31ST BIT OF FULL PRODUCT C K = FHI/B15 K = K - DMOD(K,1.0D0) C C ASSEMBLE ALL THE PARTS AND PRESUBTRACT P C THE PARENTHESES ARE ESSENTIAL C IX = (((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16) + K C C ADD P BACK IN IF NECESSARY C IF (IX.LT.0.0D0) IX = IX + P C C MULTIPLY BY 1/(2**31-1) C DRAND = IX*4.656612875D-10 RETURN C C * * * END OF FUNCTION DRAND * * * C END SUBROUTINE GINVSE(A, MA, NA, AIN, NAIN, M, N, TOL, V, NV, * FAIL, NOUT) C C INVERSE.OF.A = V * INVERSE.OF.S * TRANSPOSE.OF.U C VIA SINGULAR VALUE DECOMPOSITION USING ALGORITHMS 1 AND 2 OF C NASH J C, COMPACT NUMERICAL METHODS FOR COMPUTERS, 1979, C ADAM HILGER, BRISTOL OR HALSTED PRESS, N.Y. C C C CALLING SEQUENCE C C CALL GINVSE (A,MA,NA,AIN,NAIN,M,N,TOL,V,NV,FAIL,NOUT) C C A = MATRIX OF SIZE M BY N TO BE 'INVERTED' C A IS DESTROYED BY THIS ROUTINE AND BECOMES THE C U MATRIX OF THE SINGULAR VALUE DECOMPOSITION C MA = ROW DIMENSION OF ARRAY A C UNCHANGED BY THIS SUB-PROGRAM C NA = COLUMN DIMENSION OF ARRAY A C UNCHANGED BY THIS SUB-PROGRAM C AIN = COMPUTED 'INVERSE' (N ROWS BY M COLS.) C NAIN = ROW DIMENSION OF ARRAY AIN C ( COLUMN DIMENSION ASSUMED TO BE AT LEAST M ) C UNCHANGED BY THIS SUB-PROGRAM C M = NUMBER OF ROWS IN MATRIX A AND C THE NUMBER OF COLUMNS IN MATRIX AIN C N = NUMBER OF COLUMNS IN MATRIX A AND C THE NUMBER OF ROWS IN MATRIX AIN C UNCHANGED BY THIS SUB-PROGRAM C TOL = MACHINE PRECISION FOR COMPUTING ENVIRONMENT C UNCHANGED BY THIS SUB-PROGRAM C V = WORK MATRIX WHICH BECOMES THE V MATRIX (N BY N) C OF THE SINGULAR VALUE DECOMPOSITION C NV = DIMENSIONS OF V (BOTH ROWS AND COLUMNS) C UNCHANGED BY THIS SUB-PROGRAM C FAIL = ERROR FLAG - .TRUE. IMPLIES FAILURE OF GINVSE C NOUT = OUTPUT CHANNEL NUMBER C UNCHANGED BY THIS SUB-PROGRAM C C NOTE. THIS ROUTINE IS NOT DESIGNED TO PRODUCE A 'GOOD' C GENERALIZED INVERSE. IT IS MEANT ONLY TO FURNISH TEST C DATA FOR THE ROUTINES GMATX, ZIELKE, AND PTST C INTEGER MA, NA, NAIN, ICOUNT, N, M, NM1, I, J, JP1, K, * SWEEP, NV, NOUT, LIMIT LOGICAL FAIL REAL TOL, P, Q, R, VV, C, S, TEMP REAL A(MA,NA), AIN(NAIN,M), V(NV,NV) C C INITIALLY SET FAIL=.FALSE. TO IMPLY CORRECT EXECUTION C FAIL = .FALSE. C C TESTS ON VALIDITY OF DIMENSIONS C IF (M.LT.2 .OR. MA.LT.2 .OR. N.LT.2 .OR. NA.LT.2) GO TO * 230 IF (N.GT.NA) GO TO 210 IF (M.GT.MA) GO TO 220 C C N MUST NOT EXCEED M, OTHERWISE GINVSE WILL FAIL C IF (N.GT.M) GO TO 200 C C SWEEP COUNTER INITIALIZED TO ZERO C SWEEP = 0 C C SET SWEEP LIMIT C LIMIT = MAX0(N,6) C C MAX0(N,6) WAS CHOSEN FROM EXPERIENCE C C C V(I,J) INITIALLY N BY N IDENTITY C DO 20 I=1,N DO 10 J=1,N V(I,J) = 0.0 10 CONTINUE V(I,I) = 1.0 20 CONTINUE C C INITIALIZE ROTATION COUNTER (COUNTS DOWN TO 0) C 30 ICOUNT = N*(N-1)/2 C C COUNT SWEEP C SWEEP = SWEEP + 1 NM1 = N - 1 DO 110 J=1,NM1 JP1 = J + 1 DO 100 K=JP1,N P = 0. Q = 0. R = 0. DO 40 I=1,M C C TEST FOR AND AVOID UNDERFLOW C NOT NEEDED FOR MACHINES WHICH UNDERFLOW TO ZERO WITHOUT C ERROR MESSAGE C IF (ABS(A(I,J)).GE.TOL) Q = Q + A(I,J)*A(I,J) IF (ABS(A(I,K)).GE.TOL) R = R + A(I,K)*A(I,K) IF (ABS(A(I,J)/TOL)*ABS(A(I,K)/TOL).GE.1.0) P = P + * A(I,J)*A(I,K) 40 CONTINUE IF (Q.GE.R) GO TO 50 C = 0. S = 1. GO TO 60 50 IF (ABS(Q).LT.TOL .AND. ABS(R).LT.TOL) GO TO 90 IF (R.EQ.0.0) GO TO 90 IF ((P/Q)*(P/R).LT.TOL) GO TO 90 C C CALCULATE THE SINE AND COSINE OF THE ANGLE OF ROTATION C Q = Q - R VV = SQRT(4.*P*P+Q*Q) C = SQRT((VV+Q)/(2.*VV)) S = P/(VV*C) C C APPLY THE ROTATION TO A C 60 DO 70 I=1,M R = A(I,J) A(I,J) = R*C + A(I,K)*S A(I,K) = -R*S + A(I,K)*C 70 CONTINUE C C APPLY THE ROTATION TO V C DO 80 I=1,N R = V(I,J) V(I,J) = R*C + V(I,K)*S V(I,K) = -R*S + V(I,K)*C 80 CONTINUE GO TO 100 90 ICOUNT = ICOUNT - 1 100 CONTINUE 110 CONTINUE C C PRINT THE NUMBER OF SWEEPS AND ROTATIONS C IF (NOUT.GT.0) WRITE (NOUT,9997) SWEEP, ICOUNT C C CHECK NUMBER OF SWEEPS AND ROTATIONS (TERMINATION TEST) C IF (ICOUNT.GT.0 .AND. SWEEP.LT.LIMIT) GO TO 30 IF (SWEEP.GE.LIMIT .AND. NOUT.GT.0) WRITE (NOUT,9996) DO 160 J=1,N Q = 0. DO 120 I=1,M Q = Q + A(I,J)*A(I,J) 120 CONTINUE C C ARBITRARY RANK DECISION C IF (J.EQ.1) VV = 1.0E-3*SQRT(Q) IF (SQRT(Q).LE.VV) GO TO 140 DO 130 I=1,M A(I,J) = A(I,J)/Q 130 CONTINUE GO TO 160 140 DO 150 I=1,M A(I,J) = 0.0 150 CONTINUE 160 CONTINUE DO 190 I=1,N DO 180 J=1,M TEMP = 0.0 DO 170 K=1,N TEMP = TEMP + V(I,K)*A(J,K) 170 CONTINUE AIN(I,J) = TEMP 180 CONTINUE 190 CONTINUE RETURN C C IF AN ERROR OCCURED THE APPROPRIATE MESSAGE IS PRINTED C AND FAIL IS SET TO .TRUE. C 200 IF (NOUT.GT.0) WRITE (NOUT,9995) FAIL = .TRUE. RETURN 210 IF (NOUT.GT.0) WRITE (NOUT,9999) N, NA FAIL = .TRUE. RETURN 220 IF (NOUT.GT.0) WRITE (NOUT,9998) M, MA FAIL = .TRUE. RETURN 230 IF (NOUT.GT.0) WRITE (NOUT,9994) M, MA, N, NA FAIL = .TRUE. RETURN C C * * * * * * * * * * * * * * * * * * * C 9999 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE, * 23H FOR N IN GINVSE NA= , I5, 5(/)) 9998 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE, * 23H FOR M IN GINVSE MA= , I5, 5(/)) 9997 FORMAT (8H SWEEP =, I4, 2H , I4, 19H JACOBI ROTATIONS P, * 8HERFORMED) 9996 FORMAT (21H SWEEP LIMIT REACHED ) 9995 FORMAT (38H0**ERROR** N GREATER THAN M IS INVALID, * 21H AND GINVSE WILL FAIL, 5(/)) 9994 FORMAT (46H0**ERROR** SOME DIMENSIONS ARE LESS THAN TWO I, * 8HN GINVSE /10X, 3HM =, I5, 5H MA =, I5, 4H N =, I5, * 5H NA =, I5, 5(/)) C C * * * END OF SUBROUTINE GINVSE * * * C END SUBROUTINE GMATX(M, N, A, MA, NA, B, NB, ALPHA, K, NOUT, X ISEED, FAIL) C C FEBRUARY 9, 1984, MAY 27, 1984, JULY 18, 1985 C C BY J. C. NASH AND R. L. C. WANG C COMPUTES AN M BY N MATRIX - A - FOR USE IN TESTING PROGRAMS C PURPORTING TO CALCULATE THE GENERALIZED INVERSE OF A MATRIX C C THE MATRIX A, DIMENSION MA BY NA BUT SIZE M BY N, IS TO C BE CALCULATED BY A SEQUENCE OF PSEUDO-RANDOM JACOBI C ROTATIONS APPLIED TO A 'DIAGONAL' MATRIX WHOSE ELEMENTS C ARE DETERMINED BY THE PARAMETERS K AND ALPHA AS FOLLOWS. C A(I,I) IS SET TO ALPHA**(1-I) FOR I=1,2,...,K AND TO ZERO C FOR I.GT.K, WHERE K IS A POSITIVE INTEGER NO LARGER THAN C THE MINIMUM VALUE OF M AND N. C C CALLING SEQUENCE C C CALL GMATX(M,N,A,MA,NA,B,NB,ALPHA,K,NOUT,ISEED,FAIL) C C M = NUMBER OF ROWS IN MATRIX A C NORMALLY UNCHANGED BY THIS SUB-PROGRAM C ALTERED DURING THE EXECUTION OF THE SUB-PROGRAM IF N.GT.M C THEN RESET BEFORE RETURN TO THE CALLING PROGRAM C M IS RESET TO MA IF M.GT.MA ON ENTRY C N = NUMBER OF COLUMNS IN MATRIX A C NORMALLY UNCHANGED BY THIS SUB-PROGRAM C ALTERED DURING THE EXECUTION OF THE SUB-PROGRAM IF N.GT.M C THEN RESET BEFORE RETURN TO THE CALLING PROGRAM C NO CHECK IS MADE TO VERIFY N.LE.NA C A = THE MATRIX TO BE CREATED IN THIS SUBROUTINE C A MUST BE DECLARED WITH FIRST DIMENSION = MA C AND SECOND DIMENSION AT LEAST N IN CALLING PROGRAM C MA = FIRST OR ROW DIMENSION OF A C MA SHOULD CORRESPOND TO THE FIRST DIMENSION C OF MATRIX A IN THE CALLING PROGRAM C UNCHANGED BY THIS SUB-PROGRAM C NA = SECOND OR COLUMN DIMENSION OF MATRIX A C SHOULD NOT BE ALTERED C B = WORKING MATRIX IN WHICH TEST MATRIX IS BUILT C NB = FIRST DIMENSION OF MATRIX B C SHOULD BE AT LEAST AS LARGE AS LARGEST DIMENSION C OF MATRIX A C NOT ALTERED BY THIS SUBROUTINE C ALPHA= FACTOR USED TO GENERATE SINGULAR VALUES OF C MATRIX A. THESE ARE GENERATED ACCORDING TO C THE FORMULA ALPHA**(1-I) FOR I=1,2,...,K C WHERE K IS THE RANK (SEE BELOW) C ALPHA SHOULD LIE IN INTERVAL (0.1, 10.0) TO C GENERATE 'REASONABLE' SINGULAR VALUES, BUT C OTHER POSITIVE VALUES ARE ACCEPTED WITHOUT CHANGE C IF ALPHA .LE. 0.0 IT WILL BE SET TO 1.0 C K = THE RANK OF THE MATRIX TO BE GENERATED C MUST BE POSITIVE AND NO LARGER THAN THE C MINIMUM OF M AND N. IF K .GT. MIN(M,N) IT C WILL BE SET TO MIN(M,N) AND A WARNING C MESSAGE PRINTED C OTHERWISE K IS UNALTERED BY THIS SUB-PROGRAM C NOUT = OUTPUT CHANNEL NUMBER C SETTING NOUT=0 WILL SUPPRESS OUTPUT C UNCHANGED BY THIS SUB-PROGRAM C ISEED= SEED FOR THE PSEUDO-RANDOM NUMBER GENERATOR C ONLY POSITIVE VALUES OF ISEED ARE ALLOWED. C IF ISEED .LE. 0 IT WILL BE SET TO 1 C ISEED IS CONVERTED TO A DOUBLE PRECISION C VARIABLE DSEED FOR USE IN DRAND, THE GENERATOR C SUBPROGRAM C DSEED IS ALTERED BY DRAND, AND DSEED MUST BE C SUPPLIED IN EVERY CALL TO DRAND C REFERENCE -- C SCHRAGE L. C A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE C VOL. 5, NO. 2, JUNE 1979, PAGES 132-138 C FAIL = ERROR FLAG C FAIL = .TRUE. IF GMATX HAS FAILED TO CREATE MATRIX A C FAIL = .FALSE. IF GMATX HAS CREATED MATRIX A C MOST FAILURES ARE ASSOCIATED WITH THE USE OF C INCORRECT ARGUMENTS FOR THIS SUBROUTINE C C C NOTE -- RANDOM NUMBER GENERATOR DECLARED DOUBLE PRECISION C DOUBLE PRECISION DRAND, DSEED, DBLE INTEGER I, IM, IMAX, ISEED, J, J1, K, L, M, MA, MN, * M1, N, NA, NB, NOUT LOGICAL FAIL REAL A(MA,NA), B(NB,NB) REAL ALPHA, C, S, T C C REAL FLOAT C C ABOVE STATEMENT MAY BE NEEDED BY SOME COMPILERS C C C SET FLAG INITIALLY TO .FALSE. (PROGRAM O.K.) FAIL = .FALSE. C C ISEED IS REQUIRED TO BE POSITIVE C IF (ISEED.GT.0) GO TO 10 IF (NOUT.GT.0) WRITE (NOUT,9999) ISEED ISEED = 1 C C TESTS ON M,N AND MA C C TEST FOR FIRST DIMENSION TOO SMALL 10 IF (MA .LT. 1) GO TO 170 C TEST FOR INVALID ROW OR COLUMN SIZE IF (M .LT. 1 .OR. N .LT. 1) GO TO 170 C TEST FOR VALID ROW SIZE IF (M.LE.MA) GO TO 20 IF (NOUT.GT.0) WRITE (NOUT,9998) M, MA M = MA 20 IF (NOUT.GT.0) WRITE (NOUT,9993) IF (NOUT.GT.0) WRITE (NOUT,9992) M, N, ISEED C C DETERMINE IF A OR A-TRANSPOSE IS TO BE GENERATED BY THE C JACOBI ROTATIONS C IMAX = 1 IF A IS GENERATED, C = 0 IF A-TRANSPOSE IS GENERATED IMAX = 1 MN = N IF ( M .GE. N ) GO TO 1 IMAX = 0 IM = M M = N N = IM MN = M C 1 DO 40 I=1,M DO 30 J=1,N B(I,J) = 0.0 30 CONTINUE 40 CONTINUE C C REPLACE ( ALPHA .LE. ZERO ) BY 1.0 C IF (ALPHA.GT.0.0) GO TO 50 IF (NOUT.NE.0) WRITE (NOUT,9997) ALPHA ALPHA = 1.0 C C C TEST K AND LOOP, K BEING THE RANK C 50 IF (K.LT.1) GO TO 140 IF (K.LE.MN) GO TO 60 IF (NOUT.NE.0) WRITE (NOUT,9996) K, MN K = MN IF (NOUT.GT.0) WRITE (NOUT,9991) ALPHA, K 60 DO 70 I=1,K B(I,I) = ALPHA**(1-I) 70 CONTINUE M1 = M - 1 C C SET RANDOM NUMBER GENERATOR SEED C DSEED = DBLE(FLOAT(ISEED)) C C PERFORM JACOBI ROTATIONS (ONE SWEEP ONLY) C JPASS = 1 C C MORE SWEEPS CAN BE PERFORMED AS FOLLOWS C C JPASS = 3 C FOR THREE SWEEPS DO 120 IPASS=1,JPASS C IF (M1.LT.1) GO TO 115 C AVOID LOOP WHEN TRIVIAL 1 BY 1 MATRIX TO BE GENERATED DO 110 I=1,M1 J1 = I + 1 DO 100 J=J1,M C C GET SINE OF ANGLE OF ROTATION FROM PSEUDO RANDOM C NUMBER GENERATOR C NOTE THAT THIS IS ASSIGNED TO A SINGLE PRECISION RESULT C S = DRAND(DSEED) C C COMPUTE COSINE C C = SQRT(1.0-S*S) C C ROWS C DO 80 L=1,N T = B(I,L) B(I,L) = T*C + B(J,L)*S B(J,L) = -T*S + B(J,L)*C 80 CONTINUE C C COLUMN TRANSFORMATIONS - THESE ARE OMITTED IF C THE INDICES ARE TOO LARGE (THAT IS, IF J.GT.N) C SINCE MAX(I)=J-1, THERE IS NO TEST ON I C 85 IF (J .GT. N) GO TO 100 C DO 90 L=1,M T = B(L,I) B(L,I) = T*C + B(L,J)*S B(L,J) = -T*S + B(L,J)*C 90 CONTINUE 100 CONTINUE 110 CONTINUE C NEXT STATEMENT IS NOT THE END OF A LOOP C BUT IS THE CONTINUATION POINT AFTER LOOPS WHEN C A TRIVIAL 1 BY 1 MATRIX IS TO BE GENERATED 115 CONTINUE 120 CONTINUE C C RESET ROW AND COLUMN SIZES C IF (IMAX.EQ.1) GO TO 124 IM = M M = N N = IM C C MOVE B OR B-TRANSPOSE TO A. C 124 DO 126 I = 1,M DO 126 J=1,N IF (IMAX .EQ. 0) A(I,J) = B(J,I) IF (IMAX .EQ. 1) A(I,J) = B(I,J) 126 CONTINUE C C OUTPUT THE MATRIX WHICH HAS BEEN GENERATED C IF (NOUT.LE.0) RETURN DO 130 I=1,M WRITE (NOUT,9990) I WRITE (NOUT,9989) (A(I,J),J=1,N) 130 CONTINUE RETURN 140 IF (NOUT.GT.0) WRITE (NOUT,9995) K FAIL = .TRUE. RETURN 170 FAIL = .TRUE. IF (NOUT.GT.0) WRITE (NOUT,9994) MA, M, N RETURN C C * * * * * * * * * * * * * * * * C 9999 FORMAT (12H0**WARNING**, I10, 24H WAS AN IMPROPER CHOICE , * 40HFOR ISEED, WHICH WAS SET TO ONE IN GMATX) 9998 FORMAT (13H0**WARNING** , I5, 24H EXCEEDS ALLOWABLE VALUE, * 27H FOR M, AND WILL BE SET TO , I5, 9H IN GMATX) 9997 FORMAT (12H0**WARNING**, F10.5, 22H IS NOT ALLOWABLE FOR , * 37HALPHA AND WILL BE SET TO ONE IN GMATX) 9996 FORMAT (12H0**WARNING**, I5, 25H EXCEEDS ALLOWABLE VALUE , * 36HFOR K (THE RANK). K WILL BE SET TO , I5) 9995 FORMAT (10H0**ERROR**, I5, 27H IS AN INVALID VALUE FOR K , * 26H,THIS FORBIDS USE OF GMATX, 5(/)) 9994 FORMAT (46H0**ERROR** DIMENSIONS OF A CANNOT BE LESS THAN, * 4H 1, , 5HMA = , I5, 5H M = , I5, 5H N = , I5, 5(/)) 9993 FORMAT (42H0GENERATING MATRIX BY PSEUDO-RANDOM JACOBI, * 10H ROTATIONS) 9992 FORMAT (7H ORDER=, I5, 4H BY , I4/20H PSEUDO-RANDOM NUM, * 9HBER SEED=, I10) 9991 FORMAT (30H INITIAL DIAGONAL ELEMENTS ARE, F10.5, * 20H**(I-1) FOR I=1 TO , I3) 9990 FORMAT (4H ROW, I4) 9989 FORMAT (1H , 5(1PE16.8) ) C C * * * END OF SUBROUTINE GMATX * * * C END SUBROUTINE PTST(M, N, A, MA, X, NX, C, NC, TA, TM, NOUT, PTST 10 * FAIL) PTST 20 C PTST 30 C PTST 40 C THIS SUBROUTINE IS DESIGNED TO TEST A PROPOSED GENERALIZED PTST 50 C INVERSE OF A MATRIX LABELLED - A - WHICH IS M BY N . PTST 60 C - X - IS THE N BY M MATRIX CONTAINING THE SUPPOSED INVERSE. PTST 70 C PTST 80 C PROGRAM ORIGINALLY FOR SQUARE, SYMMETRIC MATRICES WRITTEN BY PTST 90 C RICHARD L. C. WANG, 1977 PTST 100 C MODIFIED AND ADAPTED BY J. C. NASH 1979,1982 PTST 110 C PTST 120 C PTST 130 C NOTE THAT A AND X MAY BE INTERCHANGED -- PTST 140 C TESTER DOES NOT OBJECT PTST 150 C PTST 160 C AVERAGE AND MAXIMUM ABSOLUTE DEVIATIONS ARE CALCULATED FOR PTST 170 C THE MATRICES -- PTST 180 C PTST 190 C AXA - A ... PENROSE CONDITION 1 PTST 200 C XAX - X ... PENROSE CONDITION 2 PTST 210 C (AX)T - AX ... PENROSE CONDITION 3 PTST 220 C (XA)T - XA ... PENROSE CONDITION 4 PTST 230 C PTST 240 C PTST 250 C CALLING SEQUENCE -- PTST 260 C PTST 270 C CALL PTST(M,N,A,MA,X,NX,C,NC,TA,TM,NOUT,FAIL) PTST 280 C PTST 290 C WHERE PTST 300 C PTST 310 C M = NUMBER OF ROWS IN THE 'ORIGINAL' MATRIX A PTST 320 C = NUMBER OF COLUMNS IN PURPORTED INVERSE PTST 330 C UNCHANGED BY THIS SUB-PROGRAM PTST 340 C N = NUMBER OF COLUMNS IN 'ORIGINAL' MATRIX A PTST 350 C = NUMBER OF ROWS IN PURPORTED INVERSE PTST 360 C UNCHANGED BY THIS SUB-PROGRAM PTST 370 C A = 'ORIGINAL' MATRIX OF WHICH A GENERALIZED PTST 380 C INVERSE HAS SUPPOSEDLY BEEN COMPUTED PTST 390 C UNCHANGED BY THIS SUB-PROGRAM PTST 400 C MA = FIRST OR ROW DIMENSION OF A PTST 410 C UNCHANGED BY THIS SUB-PROGRAM PTST 420 C X = PURPORTED GENERALIZED INVERSE OF MATRIX A PTST 430 C UNCHANGED BY THIS SUB-PROGRAM PTST 440 C NX = FIRST OR ROW DIMENSION OF X PTST 450 C UNCHANGED BY THIS SUB-PROGRAM PTST 460 C C = DOUBLE PRECISION WORKING ARRAY PTST 470 C SHOULD BE DIMENSIONED AT LEAST MA BY MA PTST 480 C NC = FIRST OR ROW DIMENSION OF C PTST 490 C UNCHANGED BY THIS SUB-PROGRAM PTST 500 C SHOULD BE AT LEAST AS LARGE AS MA PTST 510 C TA = VECTOR (1 DIMENSIONAL ARRAY) OF 4 ELEMENTS TO PTST 520 C STORE AVERAGE ABSOLUTE DEVIATIONS FROM EACH OF PTST 530 C THE FOUR PENROSE CONDITIONS PTST 540 C TM = VECTOR (1 DIMENSIONAL ARRAY) OF 4 ELEMENTS TO PTST 550 C STORE MAXIMUM ABSOLUTE DEVIATIONS FROM EACH OF PTST 560 C THE FOUR PENROSE CONDITIONS PTST 570 C NOUT = OUTPUT CHANNEL NUMBER PTST 580 C SETTING NOUT=0 SUPPRESSES OUTPUT PTST 590 C UNCHANGED BY THIS SUB-PROGRAM PTST 600 C FAIL = FAILURE FLAG SET .TRUE. FOR FAILURE PTST 610 C .FALSE. OTHERWISE PTST 620 C PTST 630 C THIS ROUTINE USES DOUBLE PRECISION ACCUMULATION OF PTST 640 C INNER PRODUCTS TO LIMIT ROUNDING ERROR PTST 650 C PTST 660 C REFERENCE -- WILKINSON J.H. PTST 670 C THE ALGEBRAIC EIGENVALUE PROBLEM PTST 680 C OXFORD, 1965 PTST 690 C PTST 700 C PTST 710 C CALLED FUNCTION SUBPROGRAM -- ANORM PTST 720 C PTST 730 C REAL FUNCTION ANORM(M, N, A, MA) PTST 740 C PTST 750 C ANORM COMPUTES THE SQUARE NORM OF MATRIX A PTST 760 C NORM=SQRT(SUM(A(I,J)**2), FOR I=1,M, J=1,N) PTST 770 C THIS NORM IS USED FOR SIMPLICITY PTST 780 C --- OTHER NORMS ARE ACCEPTABLE PTST 790 C PTST 800 C PTST 810 DOUBLE PRECISION DABS, DBLE, S, V PTST 820 INTEGER I, J, J1, L, M, MA, N, NC, NOUT, NX, N1, M1 PTST 830 LOGICAL FAIL PTST 840 DOUBLE PRECISION C(NC,NC) PTST 850 REAL A(MA,N), X(NX,M), TA(4), TM(4) PTST 860 REAL AMA, ANORM, ANX, T1, T2 PTST 870 C PTST 880 C DIMENSIONS MUST BE AT LEAST ONE PTST 890 C NO TESTS FOR DIMENSIONS SUFFICIENTLY LARGE ARE MADE PTST 900 C PTST 910 C INITIALIZE FAILURE FLAG TO INDICATE SUCCESSFUL OPERATION PTST 920 FAIL = .FALSE. PTST 930 IF (MA.GT.0 .AND. N.GT.0 .AND. NX.GT.0 .AND. M.GT.0 .AND. PTST 940 * NC.GT.0) GO TO 20 PTST 950 IF (NOUT.EQ.0) GO TO 10 PTST 960 WRITE (NOUT,9999) PTST 970 WRITE (NOUT,9998) MA, N, NX, M, NC PTST 980 10 FAIL = .TRUE. PTST 990 RETURN PTST1000 C PTST1010 C ZERO TEST VALUES PTST1020 C PTST1030 20 DO 30 I=1,4 PTST1040 TA(I) = 0.0 PTST1050 TM(I) = 0.0 PTST1060 30 CONTINUE PTST1070 C PTST1080 C COMPUTE AX PTST1090 C PTST1100 DO 60 I=1,M PTST1110 DO 50 J=1,M PTST1120 S = 0.0 PTST1130 DO 40 L=1,N PTST1140 S = S + A(I,L)*DBLE(X(L,J)) PTST1150 40 CONTINUE PTST1160 C(I,J) = S PTST1170 50 CONTINUE PTST1180 60 CONTINUE PTST1190 C PTST1200 C COMPUTE AXA, AXA - A = (M BY N) PTST1210 C PTST1220 DO 90 I=1,M PTST1230 DO 80 J=1,N PTST1240 S = 0.0 PTST1250 DO 70 L=1,M PTST1260 S = S + C(I,L)*DBLE(A(L,J)) PTST1270 70 CONTINUE PTST1280 T1 = DABS(DBLE(A(I,J))-S) PTST1290 IF (T1.GT.TM(1)) TM(1) = T1 PTST1300 TA(1) = TA(1) + T1 PTST1310 80 CONTINUE PTST1320 90 CONTINUE PTST1330 TA(1) = TA(1)/FLOAT(M*N) PTST1340 C PTST1350 C COMPUTE XAX, XAX - X = (N BY M) PTST1360 C PTST1370 DO 120 I=1,N PTST1380 DO 110 J=1,M PTST1390 S = 0.0 PTST1400 DO 100 L=1,M PTST1410 S = S + DBLE(X(I,L))*C(L,J) PTST1420 100 CONTINUE PTST1430 T1 = DABS(DBLE(X(I,J))-S) PTST1440 IF (T1.GT.TM(2)) TM(2) = T1 PTST1450 TA(2) = TA(2) + T1 PTST1460 110 CONTINUE PTST1470 120 CONTINUE PTST1480 TA(2) = TA(2)/FLOAT(M*N) PTST1490 C PTST1500 C ASYMMETRY OF AX (M BY M) PTST1510 C PTST1520 C TEST FOR TRIVIAL CASE PTST1530 C PTST1540 C NOTE THAT NORMALIZATION UNNECESSARY WHEN M = 1 PTST1550 C SINCE AX WILL BE 1 BY 1 PTST1560 C PTST1570 IF (M .EQ. 1) GO TO 145 PTST1580 M1 = M - 1 PTST1590 DO 140 I=1,M1 PTST1600 J1 = I + 1 PTST1610 DO 130 J=J1,M PTST1620 T1 = DABS(C(I,J)-C(J,I)) PTST1630 IF (T1.GT.TM(3)) TM(3) = T1 PTST1640 TA(3) = TA(3) + T1 PTST1650 130 CONTINUE PTST1660 140 CONTINUE PTST1670 TA(3) = TA(3)/FLOAT(M*(M-1)/2) PTST1680 C PTST1690 C ASYMMETRY OF XA (N BY N) PTST1700 C PTST1710 C TEST FOR TRIVIAL CASE PTST1720 C PTST1730 C NOTE THAT NORMALIZATION UNNECESSARY WHEN N = 1 PTST1740 C SINCE XA WILL BE 1 BY 1 PTST1750 C PTST1760 C PTST1770 145 IF (N .EQ. 1) GO TO 175 PTST1780 N1 = N - 1 PTST1790 DO 170 I=1,N1 PTST1800 J1 = I + 1 PTST1810 DO 160 J=J1,N PTST1820 S = 0.0 PTST1830 V = 0.0 PTST1840 DO 150 L=1,M PTST1850 S = S + X(I,L)*DBLE(A(L,J)) PTST1860 V = V + X(J,L)*DBLE(A(L,I)) PTST1870 150 CONTINUE PTST1880 T1 = DABS(S-V) PTST1890 IF (T1.GT.TM(4)) TM(4) = T1 PTST1900 TA(4) = TA(4) + T1 PTST1910 160 CONTINUE PTST1920 170 CONTINUE PTST1930 TA(4) = TA(4)/FLOAT(N*(N-1)/2) PTST1940 175 IF (NOUT.LE.0) RETURN PTST1950 C PTST1960 C REPORT RESULTS PTST1970 C PTST1980 WRITE (NOUT,9997) PTST1990 WRITE (NOUT,9996) PTST2000 C PTST2010 C COMPUTE NORMS IN ORDER TO GAUGE RELATIVE ERROR SIZES PTST2020 C PTST2030 AMA = ANORM(M,N,A,MA) PTST2040 ANX = ANORM(N,M,X,MA) PTST2050 WRITE (NOUT,9994) AMA, ANX PTST2060 WRITE (NOUT,9995) PTST2070 T1 = TA(1)/AMA PTST2080 T2 = TM(1)/AMA PTST2090 WRITE (NOUT,9993) TA(1), TM(1), T1, T2 PTST2100 T1 = TA(2)/ANX PTST2110 T2 = TM(2)/ANX PTST2120 WRITE (NOUT,9992) TA(2), TM(2), T1, T2 PTST2130 T1 = TA(3)/(AMA*ANX) PTST2140 T2 = TM(3)/(AMA*ANX) PTST2150 WRITE (NOUT,9991) TA(3), TM(3), T1, T2 PTST2160 T1 = TA(4)/(AMA*ANX) PTST2170 T2 = TM(4)/(AMA*ANX) PTST2180 WRITE (NOUT,9990) TA(4), TM(4), T1, T2 PTST2190 RETURN PTST2200 C PTST2210 C * * * * * * * * * * * * * * * * PTST2220 C PTST2230 9999 FORMAT (46H0**ERROR** SOME DIMENSIONS ARE LESS THAN ONE I, PTST2240 * 6HN PTST) PTST2250 9998 FORMAT (11H MA=, I5, 3H N=, I5, 4H NX=, I5, 3H M=, PTST2260 * I5, 4H NC=, I5, 5(/)) PTST2270 9997 FORMAT (25H0 TEST PENROSE CONDITIONS/) PTST2280 9996 FORMAT (46H A IS THE INPUT MATRIX, X IS THE INVERSE, PTST2290 * 5H OF A) PTST2300 9995 FORMAT (1H , 35X, 36HAVERAGE DEVIATION MAXIMUM DEVIATI, PTST2310 * 2HON) PTST2320 9994 FORMAT (22H0 INPUT MATRIX NORM= , 1PE16.8/9H INVERSE, PTST2330 * 13H MATRIX NORM=, 1PE16.8) PTST2340 9993 FORMAT (29H0 AXA=TEST=A ... ACTUAL , 7X, 2(E16.8, PTST2350 * 5X)/19X, 10HNORMALIZED, 7X, 2(E16.8, 5X)) PTST2360 9992 FORMAT (29H0 XAX=TEST=X ... ACTUAL , 7X, 2(E16.8, PTST2370 * 5X)/19X, 10HNORMALIZED, 7X, 2(E16.8, 5X)) PTST2380 9991 FORMAT (29H0 (AX)T=TEST=AX ... ACTUAL , 7X, 2(E16.8, PTST2390 * 5X)/19X, 10HNORMALIZED, 7X, 2(E16.8, 5X)) PTST2400 9990 FORMAT (29H0 (XA)T=TEST=XA ... ACTUAL , 7X, 2(E16.8, PTST2410 * 5X)/19X, 10HNORMALIZED, 7X, 2(E16.8, 5X), 5(/)) PTST2420 C PTST2430 C * * * END OF SUBROUTINE PTST * * * PTST2440 C PTST2450 END PTST2460 REAL FUNCTION ANORM(M, N, A, MA) PTST2470 C PTST2480 C COMPUTE SQUARE (EUCLIDEAN) NORM OF MATRIX A PTST2490 C NORM=SQRT(SUM(A(I,J)**2), FOR I=1,M, J=1,N) PTST2500 C THIS NORM IS USED FOR SIMPLICITY - OTHERS ACCEPTABLE PTST2510 C PTST2520 C CALLING SEQUENCE PTST2530 C PTST2540 C RNORM = ANORM(M,N,A,MA) PTST2550 C PTST2560 C RNORM = REAL VARIABLE WHICH RETURNS THE VALUE OF THE NORM PTST2570 C VIA THE FUNCTION PTST2580 C M = NUMBER OF ROWS IN MATRIX A PTST2590 C N = NUMBER OF COLUMNS IN MATRIX A PTST2600 C A = SUBJECT MATRIX WHICH IS M BY N PTST2610 C MA = FIRST DIMENSION OF A PTST2620 C PTST2630 C NONE OF THESE ARGUMENTS ARE ALTERED BY THIS SUB-PROGRAM PTST2640 C PTST2650 C DOUBLE PRECISION ACCUMULATION IS USED TO LIMIT ROUNDING ERROR PTST2660 C PTST2670 INTEGER I, J, M, MA, N PTST2680 REAL A(MA,N) PTST2690 DOUBLE PRECISION DNORM, DBLE, DSQRT PTST2700 DNORM = 0.0D0 PTST2710 DO 20 I=1,M PTST2720 DO 10 J=1,N PTST2730 DNORM = DNORM + DBLE(A(I,J))**2 PTST2740 10 CONTINUE PTST2750 20 CONTINUE PTST2760 C PTST2770 C NOTE IMPLICIT CONVERSION TO SINGLE PRECISION PTST2780 C PTST2790 ANORM = DSQRT(DNORM) PTST2800 RETURN PTST2810 C PTST2820 C * * * END OF FUNCTION ANORM * * * PTST2830 C * * * END OF COMBINED SET -- PTST AND ANORM * * * PTST2840 C PTST2850 END PTST2860 SUBROUTINE ZIELKE(M, N, A, MA, NA, MOPT, ALPHA, NOUT, ZIEL 10 * FAIL) ZIEL 20 C ZIEL 30 C COMPUTES GENERALIZED INVERSE TEST MATRICES DUE TO ZIEL 40 C ZIELKE, SIGNUM NEWSLETTER, VOL. 13, #4, DEC. 78 ZIEL 50 C PAGES 10-12 ZIEL 60 C AND ZIEL 70 C ZIELKE G., SIGNUM NEWSLETTER, VOL 16, #3, SEPT. 81, ZIEL 80 C PAGES 7-8 ZIEL 90 C NOTE CORRECTIONS IN SIGNUM NEWSLETTER, VOL. 16, ZIEL 100 C #4, DEC. 81, PAGE 6 ZIEL 110 C THESE MATRICES ARE LABELLED A1, A2, A3, A4 ZIEL 120 C OR THEIR GENERALIZED INVERSES A1+, A2+, A3+, A4+ ZIEL 130 C ZIEL 140 C CALLING SEQUENCE -- ZIEL 150 C ZIEL 160 C CALL ZIELKE(M,N,A,MA,NA,MOPT,ALPHA,NOUT,FAIL) ZIEL 170 C ZIEL 180 C M = NUMBER OF ROWS IN MATRIX A PRODUCED ZIEL 190 C THIS IS SET (CHANGED) BY THIS SUBROUTINE ZIEL 200 C N = NUMBER OF COLUMNS IN MATRIX A PRODUCED ZIEL 210 C THIS IS SET (CHANGED) BY THIS SUBROUTINE ZIEL 220 C A = THE MATRIX WHICH IS GENERATED ZIEL 230 C MA = THE FIRST OR ROW DIMENSION OF A IN THE CALLING ZIEL 240 C PROGRAM ZIEL 250 C UNCHANGED BY THIS SUB-PROGRAM ZIEL 260 C NA = THE SECOND OR COLUMN DIMENSION OF A IN THE ZIEL 270 C CALLING PROGRAM ZIEL 280 C UNCHANGED BY THIS SUB-PROGRAM ZIEL 290 C MOPT = AN INTEGER USED TO SELECT THE MATRIX TO BE ZIEL 300 C GENERATED. ZIEL 310 C UNCHANGED BY THIS SUB-PROGRAM ZIEL 320 C THE FOLLOWING TABLE DETERMINES THE POSSIBLE ZIEL 330 C MATRICES WHICH MAY BE GENERATED -- ZIEL 340 C MOPT=-4 YIELDS ZIELKE/RUTISHAUSER A4+ MATRIX ZIEL 350 C (THE MOORE-PENROSE INVERSE OF A4) ZIEL 360 C MOPT=-3 YIELDS ZIELKE A3+ MATRIX ZIEL 370 C MOPT=-2 YIELDS ZIELKE A2+ MATRIX ZIEL 380 C MOPT=-1 YIELDS ZIELKE A1+ MATRIX ZIEL 390 C MOPT= 1 YIELDS ZIELKE A1 MATRIX ZIEL 400 C MOPT= 2 YIELDS ZIELKE A2 MATRIX ZIEL 410 C MOPT= 3 YIELDS ZIELKE A3 MATRIX ZIEL 420 C MOPT= 4 YIELDS ZIELKE/RUTISHAUSER A4 MATRIX ZIEL 430 C ZIEL 440 C ALPHA = PARAMETER USED IN GENERATING THE ZIELKE ZIEL 450 C MATRICES ZIEL 460 C UNCHANGED BY THIS SUB-PROGRAM ZIEL 470 C NOUT = OUTPUT CHANNEL NUMBER. ZIEL 480 C SETTING NOUT=0 SUPPRESSES OUTPUT. ZIEL 490 C UNCHANGED BY THIS SUB-PROGRAM ZIEL 500 C FAIL = FAILURE FLAG SET .TRUE. ON FAILURE, ZIEL 510 C .FALSE. OTHERWISE ZIEL 520 C ZIEL 530 C ZIEL 540 C ZIEL 550 C THE DECLARED DIMENSIONS MUST BE LARGE ENOUGH TO HOLD THE ZIEL 560 C RESULTING MATRICES. THE FOLLOWING TABLE IS USEFUL -- ZIEL 570 C MOPT MIN. MA MIN. NA ZIEL 580 C 1 5 4 ZIEL 590 C 2 5 4 ZIEL 600 C 3 5 4 ZIEL 610 C 4 7 3 ZIEL 620 C -1 4 5 ZIEL 630 C -2 4 5 ZIEL 640 C -3 4 5 ZIEL 650 C -4 3 7 ZIEL 660 C ZIEL 670 INTEGER I, J, M, MA, MOPT, N, NA, NOPT, NOUT ZIEL 680 LOGICAL FAIL ZIEL 690 REAL A(MA,NA) ZIEL 700 REAL ALPHA ZIEL 710 C ZIEL 720 C INITIALLY FAILURE FLAG INDICATES SUCCESSFUL OPERATION ZIEL 730 FAIL = .FALSE. ZIEL 740 C ZIEL 750 C USE THE VALUE OF MOPT TO DETERMINE ZIEL 760 C WHICH INPUT MATRIX IS DESIRED ZIEL 770 C ZIEL 780 NOPT = MOPT + 5 ZIEL 790 C ZIEL 800 C SPACE VERTICALLY BEFORE PRINTING ZIEL 810 C ZIEL 820 IF (NOUT.GT.0) WRITE (NOUT,9997) ZIEL 830 C ZIEL 840 C SAFETY CHECK ZIEL 850 C ZIEL 860 IF (NOPT.LT.1 .OR. NOPT.GT.9) GO TO 190 ZIEL 870 GO TO (10, 20, 30, 40, 190, 50, 60, 100, 140), NOPT ZIEL 880 C ZIEL 890 C COMPUTE ZIELKE RUTISHAUSER MATRIX A4+ ZIEL 900 C ZIEL 910 10 IF (MA.LT.3 .OR. NA.LT.7) GO TO 210 ZIEL 920 A(1,1) = -(9.*ALPHA+63.)/168. ZIEL 930 A(1,2) = -(6.*ALPHA-46.)/168. ZIEL 940 A(1,3) = -(3.*ALPHA-29.)/168. ZIEL 950 A(1,4) = -12./168. ZIEL 960 A(1,5) = (3.*ALPHA+5.)/168. ZIEL 970 A(1,6) = (6.*ALPHA+22.)/168. ZIEL 980 A(1,7) = (9.*ALPHA+39.)/168. ZIEL 990 A(2,1) = -6./168. ZIEL1000 A(2,2) = -4./168. ZIEL1010 A(2,3) = -2./168. ZIEL1020 A(2,4) = 0./168. ZIEL1030 A(2,5) = 2./168. ZIEL1040 A(2,6) = 4./168. ZIEL1050 A(2,7) = 6./168. ZIEL1060 A(3,1) = (9.*ALPHA+51.)/168. ZIEL1070 A(3,2) = (6.*ALPHA+38.)/168. ZIEL1080 A(3,3) = (3.*ALPHA+25.)/168. ZIEL1090 A(3,4) = 12./168. ZIEL1100 A(3,5) = -(3.*ALPHA+1.)/168. ZIEL1110 A(3,6) = -(6.*ALPHA+14.)/168. ZIEL1120 A(3,7) = -(9.*ALPHA+27.)/168. ZIEL1130 M = 3 ZIEL1140 N = 7 ZIEL1150 IF (NOUT.GT.0) WRITE (NOUT,9996) ZIEL1160 GO TO 170 ZIEL1170 C ZIEL1180 C COMPUTE ZIELKE MATRIX A3+ ZIEL1190 C ZIEL1200 20 IF (MA.LT.4 .OR. NA.LT.5) GO TO 210 ZIEL1210 A(1,1) = 0.5 ZIEL1220 A(1,2) = -.125 ZIEL1230 A(1,3) = -1. ZIEL1240 A(1,4) = .875 ZIEL1250 A(1,5) = -.625 ZIEL1260 A(1,6) = .375 ZIEL1270 A(2,1) = -1. ZIEL1280 A(2,2) = (2.*ALPHA+13.)/8. ZIEL1290 A(2,3) = (-8.*ALPHA-28.)/8. ZIEL1300 A(2,4) = (6.*ALPHA+17.)/8. ZIEL1310 A(2,5) = (-2.*ALPHA-3.)/8. ZIEL1320 A(2,6) = (2.*ALPHA+1.)/8. ZIEL1330 A(3,1) = 1.25 ZIEL1340 A(3,2) = -A(2,2) + .25 ZIEL1350 A(3,3) = -A(2,3) - 1.25 ZIEL1360 A(3,4) = -A(2,4) + 1. ZIEL1370 A(3,5) = -A(2,5) - .5 ZIEL1380 A(3,6) = -A(2,6) + .25 ZIEL1390 A(4,1) = -.25 ZIEL1400 A(4,2) = .375 ZIEL1410 A(4,3) = -.25 ZIEL1420 A(4,4) = .125 ZIEL1430 A(4,5) = .125 ZIEL1440 A(4,6) = -.125 ZIEL1450 A(5,1) = -.5 ZIEL1460 A(5,2) = -.25 ZIEL1470 A(5,3) = 1.5 ZIEL1480 A(5,4) = -1.25 ZIEL1490 A(5,5) = .75 ZIEL1500 A(5,6) = -.25 ZIEL1510 M = 5 ZIEL1520 N = 6 ZIEL1530 IF (NOUT.GT.0) WRITE (NOUT,9995) ZIEL1540 GO TO 170 ZIEL1550 C ZIEL1560 C COMPUTE ZIELKE MATRIX A2+ ZIEL1570 C ZIEL1580 30 IF (MA.LT.4 .OR. NA.LT.5) GO TO 210 ZIEL1590 A(1,1) = (12.*ALPHA+44.)/60. ZIEL1600 A(1,2) = 1./3. ZIEL1610 A(1,3) = (-12.*ALPHA-4.)/60. ZIEL1620 A(1,4) = (-6.*ALPHA-27.)/60. ZIEL1630 A(1,5) = (6.*ALPHA-3.)/60. ZIEL1640 A(2,1) = -A(1,1) - .2 ZIEL1650 A(2,2) = -A(1,2) ZIEL1660 A(2,3) = -A(1,3) + .2 ZIEL1670 A(2,4) = -A(1,4) + .1 ZIEL1680 A(2,5) = -A(1,5) - .1 ZIEL1690 A(3,1) = A(2,1) + 11./15. ZIEL1700 A(3,2) = 0.0 ZIEL1710 A(3,3) = A(2,3) - 1./15. ZIEL1720 A(3,4) = A(2,4) - .7 ZIEL1730 A(3,5) = A(2,5) - .3 ZIEL1740 A(4,1) = ALPHA/5. ZIEL1750 A(4,2) = 0.0 ZIEL1760 A(4,3) = -ALPHA/5. ZIEL1770 A(4,4) = -A(3,4) + .1 ZIEL1780 A(4,5) = -A(3,5) - .1 ZIEL1790 M = 4 ZIEL1800 N = 5 ZIEL1810 IF (NOUT.GT.0) WRITE (NOUT,9994) ZIEL1820 GO TO 170 ZIEL1830 C ZIEL1840 C COMPUTE ZIELKE MATRIX A1+ ZIEL1850 C ZIEL1860 40 IF (MA.LT.4 .OR. NA.LT.5) GO TO 210 ZIEL1870 A(1,1) = ALPHA/2. ZIEL1880 A(1,2) = .5 ZIEL1890 A(1,3) = ALPHA/2. ZIEL1900 A(1,4) = .5 ZIEL1910 A(1,5) = -ALPHA ZIEL1920 A(2,1) = 0.0 ZIEL1930 A(2,2) = -.25 ZIEL1940 A(2,3) = 0.0 ZIEL1950 A(2,4) = -.25 ZIEL1960 A(2,5) = .5 ZIEL1970 A(3,1) = -2.*(ALPHA+1.)/4. ZIEL1980 A(3,2) = 0.0 ZIEL1990 A(3,3) = A(3,1) ZIEL2000 A(3,4) = 0.0 ZIEL2010 A(3,5) = ALPHA ZIEL2020 A(4,1) = 0.0 ZIEL2030 A(4,2) = -.25 ZIEL2040 A(4,3) = 0.0 ZIEL2050 A(4,4) = -.25 ZIEL2060 A(4,5) = .5 ZIEL2070 M = 4 ZIEL2080 N = 5 ZIEL2090 IF (NOUT.GT.0) WRITE (NOUT,9993) ZIEL2100 GO TO 170 ZIEL2110 C ZIEL2120 C COMPUTE ZIELKE MATRIX A1 ZIEL2130 C ZIEL2140 50 IF (MA.LT.5 .OR. NA.LT.4) GO TO 210 ZIEL2150 A(1,1) = ALPHA ZIEL2160 A(1,2) = ALPHA ZIEL2170 A(1,3) = ALPHA - 1. ZIEL2180 A(1,4) = ALPHA ZIEL2190 A(2,1) = ALPHA + 1. ZIEL2200 A(2,2) = ALPHA ZIEL2210 A(2,3) = ALPHA ZIEL2220 A(2,4) = ALPHA ZIEL2230 A(3,1) = ALPHA ZIEL2240 A(3,2) = ALPHA ZIEL2250 A(3,3) = ALPHA - 1. ZIEL2260 A(3,4) = ALPHA ZIEL2270 A(4,1) = ALPHA + 1. ZIEL2280 A(4,2) = ALPHA ZIEL2290 A(4,3) = ALPHA ZIEL2300 A(4,4) = ALPHA ZIEL2310 A(5,1) = ALPHA + 1. ZIEL2320 A(5,2) = ALPHA + 1. ZIEL2330 A(5,3) = ALPHA ZIEL2340 A(5,4) = ALPHA + 1. ZIEL2350 M = 5 ZIEL2360 N = 4 ZIEL2370 IF (NOUT.GT.0) WRITE (NOUT,9992) ZIEL2380 GO TO 170 ZIEL2390 C ZIEL2400 C COMPUTE ZIELKE MATRIX A2 ZIEL2410 C ZIEL2420 60 IF (MA.LT.5 .OR. NA.LT.4) GO TO 210 ZIEL2430 A(1,1) = ALPHA + 1. ZIEL2440 A(1,2) = ALPHA ZIEL2450 A(1,3) = ALPHA ZIEL2460 A(1,4) = ALPHA + 1. ZIEL2470 DO 70 I=1,4 ZIEL2480 A(2,I) = A(1,I) + 1. ZIEL2490 70 CONTINUE ZIEL2500 DO 80 I=1,4 ZIEL2510 A(3,I) = A(2,I) + 1. ZIEL2520 80 CONTINUE ZIEL2530 A(4,1) = ALPHA + 1. ZIEL2540 A(4,2) = ALPHA + 1. ZIEL2550 A(4,3) = ALPHA ZIEL2560 A(4,4) = ALPHA + 2. ZIEL2570 DO 90 I=1,4 ZIEL2580 A(5,I) = A(4,I) - 1. ZIEL2590 90 CONTINUE ZIEL2600 M = 5 ZIEL2610 N = 4 ZIEL2620 IF (NOUT.GT.0) WRITE (NOUT,9991) ZIEL2630 GO TO 170 ZIEL2640 C ZIEL2650 C COMPUTE ZIELKE MATRIX A3 ZIEL2660 C ZIEL2670 100 IF (MA.LT.6 .OR. NA.LT.5) GO TO 210 ZIEL2680 A(1,1) = ALPHA ZIEL2690 A(1,2) = ALPHA + 1. ZIEL2700 A(1,3) = ALPHA + 2. ZIEL2710 A(1,4) = ALPHA + 3. ZIEL2720 A(1,5) = ALPHA ZIEL2730 A(2,1) = ALPHA ZIEL2740 A(2,2) = ALPHA + 2. ZIEL2750 A(2,3) = ALPHA + 3. ZIEL2760 A(2,4) = ALPHA + 5. ZIEL2770 A(2,5) = ALPHA + 1. ZIEL2780 DO 110 I=1,4 ZIEL2790 A(3,I) = ALPHA + FLOAT(I) ZIEL2800 110 CONTINUE ZIEL2810 A(3,5) = ALPHA + 2. ZIEL2820 DO 120 I=1,5 ZIEL2830 A(4,I) = A(3,I) + 1. ZIEL2840 120 CONTINUE ZIEL2850 DO 130 I=1,4 ZIEL2860 A(5,I) = A(4,I) + 1. ZIEL2870 130 CONTINUE ZIEL2880 A(5,5) = ALPHA + 5. ZIEL2890 A(6,1) = ALPHA + 5. ZIEL2900 A(6,2) = ALPHA + 5. ZIEL2910 A(6,3) = ALPHA + 6. ZIEL2920 A(6,4) = ALPHA + 6. ZIEL2930 A(6,5) = ALPHA + 7. ZIEL2940 M = 6 ZIEL2950 N = 5 ZIEL2960 IF (NOUT.GT.0) WRITE (NOUT,9990) ZIEL2970 GO TO 170 ZIEL2980 C ZIEL2990 C COMPUTE ZIELKE/RUTISHAUSER MATRIX A4 ZIEL3000 C ZIEL3010 140 IF (MA.LT.7 .OR. NA.LT.3) GO TO 210 ZIEL3020 DO 160 I=1,7 ZIEL3030 DO 150 J=1,3 ZIEL3040 A(I,J) = ALPHA + FLOAT(I+J-1) ZIEL3050 150 CONTINUE ZIEL3060 160 CONTINUE ZIEL3070 M = 7 ZIEL3080 N = 3 ZIEL3090 IF (NOUT.GT.0) WRITE (NOUT,9989) ZIEL3100 C ZIEL3110 C OUTPUT MATRIX CHOSEN ZIEL3120 C ZIEL3130 C ZIEL3140 C NORMAL EXIT FROM SUBROUTINE ZIEL3150 C ZIEL3160 170 IF (NOUT.EQ.0) RETURN ZIEL3170 DO 180 I=1,M ZIEL3180 IF (NOUT.GT.0) WRITE (NOUT,9988) I ZIEL3190 IF (NOUT.GT.0) WRITE (NOUT,9987) (A(I,J),J=1,N) ZIEL3200 180 CONTINUE ZIEL3210 RETURN ZIEL3220 C ZIEL3230 C ERROR EXIT -- WRITE APPROPRIATE MESSAGE, SET FAIL TO .TRUE. ZIEL3240 C ZIEL3250 190 IF (NOUT.EQ.0) GO TO 200 ZIEL3260 WRITE (NOUT,9999) MOPT ZIEL3270 200 FAIL = .TRUE. ZIEL3280 RETURN ZIEL3290 210 IF (NOUT.EQ.0) GO TO 220 ZIEL3300 WRITE (NOUT,9998) MA, NA, MOPT ZIEL3310 220 FAIL = .TRUE. ZIEL3320 RETURN ZIEL3330 C ZIEL3340 C * * * * * * * * * * * * * * * ZIEL3350 C ZIEL3360 9999 FORMAT (17H0**ERROR** OPTION, I5, 20H IS NOT AVAILABLE FO, ZIEL3370 * 8HR ZIELKE, 5(/)) ZIEL3380 9998 FORMAT (16H0**ERROR** MA = , I3, 10H AND NA = , I3, ZIEL3390 * 16H FORBIDS OPTION , I2, 10H IN ZIELKE, 5(/)) ZIEL3400 9997 FORMAT (1H0) ZIEL3410 9996 FORMAT (30H ZIELKE/RUTISHAUSER MATRIX A4+) ZIEL3420 9995 FORMAT (18H ZIELKE MATRIX A3+) ZIEL3430 9994 FORMAT (18H ZIELKE MATRIX A2+) ZIEL3440 9993 FORMAT (18H ZIELKE MATRIX A1+) ZIEL3450 9992 FORMAT (17H ZIELKE MATRIX A1) ZIEL3460 9991 FORMAT (17H ZIELKE MATRIX A2) ZIEL3470 9990 FORMAT (17H ZIELKE MATRIX A3) ZIEL3480 9989 FORMAT (29H ZIELKE/RUTISHAUSER MATRIX A4) ZIEL3490 9988 FORMAT (4H ROW, I4) ZIEL3500 9987 FORMAT (1H , 5(1PE16.8) ) ZIEL3510 C ZIEL3520 C * * * END OF SUBROUTINE ZIELKE * * * ZIEL3530 C ZIEL3540 END ZIEL3550 YES, DO TRY CALLING SEQUENCE TEST YES, DO TRY PTST FIXED TESTS TESTING INVALID OPTION A9+ 0 0 0 -9 0 2.0 TESTING ISEED.LT.0 IN GMATX 3 3 2 0 -1 0.2 TESTING M.LT.1 IN GMATX -1 1 1 0 34 0.1 TESTING ALPHA.LT.0.0 IN GMATX 4 4 4 0 56 -1.8 TESTING K.LT.0 IN GMATX 5 5 -1 0 98 .09 TESTING K.GT.MIN0(M,N) IN GMATX 5 3 8 0 1 1.0 TESTING N.GT.M IN GMATX, GINVSE NOT CALLED 3 4 1 0 0 1.0 TEST NO ERROR DATA SET WHICH GENERATES IDENTITY MATRIX 4 4 4 0 453 1.0 TEST 1 BY 5 ROW VECTOR 1 5 1 0 12345 0.5 TEST 5 BY 1 COLUMN VECTOR 5 1 1 0 12345 0.5 TEST 1 BY 1 TRIVIAL CASE 1 1 1 0 12345 0.5 TEST 10 BY 5 MATRIX WITH POWER OF 2 SINGULAR VALUES 10 5 5 0 12345 0.5 TEST 10 BY 5 GENERATED MATRIX WITH RANK 2 10 5 2 0 54321 0.5 TEST 5 BY 10 GENERATED MATRIX -- POWER OF 0.1, RANK 3 5 10 3 0 34521 0.1 TEST ZIELKE A1 5 5 5 1 0 0.1 TEST ZIELKE A2 5 5 5 2 0 0.1 TEST ZIELKE A3 5 5 5 3 0 0.1 TEST ZIELKE / RUTISHAUSER A4 5 5 5 4 0 0.1 TEST ZIELKE INVERSE A1 5 5 5 -1 0 0.1 TEST ZIELKE INVERSE A2 5 5 5 -2 0 0.1 TEST ZIELKE INVERSE A3 5 5 5 -3 0 0.1 TEST ZIELKE / RUTISHAUSER A4 INVERSE 5 5 5 -4 0 0.1 END OF TEST DATA SET - BLANK CARD IMAGE NO TEST OF CALLING SEQUENCES NO TESTS OF PTST WITH SMALL MATRICES SAMPLE RUN FOR ZIELKE MATRIX A1 1 9.0 SAMPLE DATA FOR GENERATING MATRIX WITH GMATX 4 4 4 123453 9.0 END OF SAMPLE DATA SET - BLANK CARD IMAGE