root / src / Ginvse.f90 @ 2
Historique | Voir | Annoter | Télécharger (8,53 ko)
| 1 |
|
|---|---|
| 2 |
SUBROUTINE GINVSE(A, MA, NA, AIN, NAIN, M, N, TOL, V, NV, FAIL, NOUT) |
| 3 |
|
| 4 |
Use Io_module |
| 5 |
! INVERSE.OF.A = V * INVERSE.OF.S * TRANSPOSE.OF.U |
| 6 |
! VIA SINGULAR VALUE DECOMPOSITION USING ALGORITHMS 1 AND 2 OF |
| 7 |
! NASH J C, COMPACT NUMERICAL METHODS FOR COMPUTERS, 1979, |
| 8 |
! ADAM HILGER, BRISTOL OR HALSTED PRESS, N.Y. |
| 9 |
! |
| 10 |
! |
| 11 |
! CALLING SEQUENCE |
| 12 |
! |
| 13 |
! CALL GINVSE (A,MA,NA,AIN,NAIN,M,N,TOL,V,NV,FAIL,NOUT) |
| 14 |
! |
| 15 |
! A = MATRIX OF SIZE M BY N TO BE 'INVERTED' |
| 16 |
! A IS DESTROYED BY THIS ROUTINE AND BECOMES THE |
| 17 |
! U MATRIX OF THE SINGULAR VALUE DECOMPOSITION |
| 18 |
! MA = ROW DIMENSION OF ARRAY A |
| 19 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 20 |
! NA = COLUMN DIMENSION OF ARRAY A |
| 21 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 22 |
! AIN = COMPUTED 'INVERSE' (N ROWS BY M COLS.) |
| 23 |
! NAIN = ROW DIMENSION OF ARRAY AIN |
| 24 |
! ( COLUMN DIMENSION ASSUMED TO BE AT LEAST M ) |
| 25 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 26 |
! M = NUMBER OF ROWS IN MATRIX A AND |
| 27 |
! THE NUMBER OF COLUMNS IN MATRIX AIN |
| 28 |
! N = NUMBER OF COLUMNS IN MATRIX A AND |
| 29 |
! THE NUMBER OF ROWS IN MATRIX AIN |
| 30 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 31 |
! TOL = MACHINE PRECISION FOR COMPUTING ENVIRONMENT |
| 32 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 33 |
! V = WORK MATRIX WHICH BECOMES THE V MATRIX (N BY N) |
| 34 |
! OF THE SINGULAR VALUE DECOMPOSITION |
| 35 |
! NV = DIMENSIONS OF V (BOTH ROWS AND COLUMNS) |
| 36 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 37 |
! FAIL = ERROR FLAG - .TRUE. IMPLIES FAILURE OF GINVSE |
| 38 |
! NOUT = OUTPUT CHANNEL NUMBER |
| 39 |
! UNCHANGED BY THIS SUB-PROGRAM |
| 40 |
! |
| 41 |
! NOTE. THIS ROUTINE IS NOT DESIGNED TO PRODUCE A 'GOOD' |
| 42 |
! GENERALIZED INVERSE. IT IS MEANT ONLY TO FURNISH TEST |
| 43 |
! DATA FOR THE ROUTINES GMATX, ZIELKE, AND PTST |
| 44 |
! |
| 45 |
INTEGER(KINT), INTENT(IN) :: MA, NA, NAIN,N, M,NV, NOUT |
| 46 |
INTEGER(KINT) :: ICOUNT, NM1, I, J, JP1, K,SWEEP, LIMIT |
| 47 |
LOGICAL, INTENT(INOUT) :: FAIL |
| 48 |
REAL(KREAL), INTENT(IN) :: TOL |
| 49 |
REAL(KREAL) :: P, Q, R, VV, C, S, TEMP |
| 50 |
REAL(KREAL), INTENT(INOUT) :: V(NV,NV) |
| 51 |
REAL(KREAL), INTENT(INOUT) :: AIN(NAIN,M) |
| 52 |
REAL(KREAL), INTENT(INOUT) :: A(MA,NA) |
| 53 |
|
| 54 |
LOGICAL :: debug |
| 55 |
INTERFACE |
| 56 |
function valid(string) result (isValid) |
| 57 |
CHARACTER(*), intent(in) :: string |
| 58 |
logical :: isValid |
| 59 |
END function VALID |
| 60 |
|
| 61 |
END INTERFACE |
| 62 |
|
| 63 |
debug=valid('Ginvse')
|
| 64 |
if (debug) WRITE(*,*) '=======Entering Ginvse=======' |
| 65 |
|
| 66 |
!Print *, 'A=' |
| 67 |
DO J=1,MA |
| 68 |
!WRITE(*,'(50(1X,F10.2))') AIN(:,J) |
| 69 |
!Print *, A(:,J) |
| 70 |
END DO |
| 71 |
!Print *, 'End of A' |
| 72 |
! |
| 73 |
! INITIALLY SET FAIL=.FALSE. TO IMPLY CORRECT EXECUTION |
| 74 |
! |
| 75 |
FAIL = .FALSE. |
| 76 |
! |
| 77 |
! TESTS ON VALIDITY OF DIMENSIONS |
| 78 |
! |
| 79 |
IF (M.LT.2 .OR. MA.LT.2 .OR. N.LT.2 .OR. NA.LT.2) GO TO 230 |
| 80 |
IF (N.GT.NA) GO TO 210 |
| 81 |
IF (M.GT.MA) GO TO 220 |
| 82 |
! |
| 83 |
! N MUST NOT EXCEED M, OTHERWISE GINVSE WILL FAIL |
| 84 |
! |
| 85 |
IF (N.GT.M) GO TO 200 |
| 86 |
! |
| 87 |
! SWEEP COUNTER INITIALIZED TO ZERO |
| 88 |
! |
| 89 |
SWEEP = 0 |
| 90 |
! |
| 91 |
! SET SWEEP LIMIT |
| 92 |
! |
| 93 |
LIMIT = MAX0(N,6) |
| 94 |
! |
| 95 |
! MAX0(N,6) WAS CHOSEN FROM EXPERIENCE |
| 96 |
! |
| 97 |
! |
| 98 |
! V(I,J) INITIALLY N BY N IDENTITY |
| 99 |
! |
| 100 |
DO 20 I=1,N |
| 101 |
DO 10 J=1,N |
| 102 |
V(I,J) = 0.0 |
| 103 |
10 CONTINUE |
| 104 |
V(I,I) = 1.0 |
| 105 |
20 CONTINUE |
| 106 |
! |
| 107 |
! INITIALIZE ROTATION COUNTER (COUNTS DOWN TO 0) |
| 108 |
! |
| 109 |
30 ICOUNT = N*(N-1)/2 |
| 110 |
! |
| 111 |
! COUNT SWEEP |
| 112 |
! |
| 113 |
SWEEP = SWEEP + 1 |
| 114 |
NM1 = N - 1 |
| 115 |
DO 110 J=1,NM1 |
| 116 |
JP1 = J + 1 |
| 117 |
DO 100 K=JP1,N |
| 118 |
P = 0. |
| 119 |
Q = 0. |
| 120 |
R = 0. |
| 121 |
DO 40 I=1,M |
| 122 |
! |
| 123 |
! TEST FOR AND AVOID UNDERFLOW |
| 124 |
! NOT NEEDED FOR MACHINES WHICH UNDERFLOW TO ZERO WITHOUT |
| 125 |
! ERROR MESSAGE |
| 126 |
! |
| 127 |
IF (ABS(A(I,J)).GE.TOL) Q = Q + A(I,J)*A(I,J) |
| 128 |
!Print *, '1 Q=',Q |
| 129 |
!Print *, 'A(',I,J,')=',A(I,J)
|
| 130 |
IF (ABS(A(I,K)).GE.TOL) R = R + A(I,K)*A(I,K) |
| 131 |
IF (ABS(A(I,J)/TOL)*ABS(A(I,K)/TOL).GE.1.0) P = P + A(I,J)*A(I,K) |
| 132 |
40 CONTINUE |
| 133 |
!Print *, '2 Q=',Q |
| 134 |
IF (Q.GE.R) GO TO 50 |
| 135 |
C = 0. |
| 136 |
S = 1. |
| 137 |
GO TO 60 |
| 138 |
50 IF (ABS(Q).LT.TOL .AND. ABS(R).LT.TOL) GO TO 90 |
| 139 |
IF (R.EQ.0.0) GO TO 90 |
| 140 |
IF ((P/Q)*(P/R).LT.TOL) GO TO 90 |
| 141 |
! |
| 142 |
! CALCULATE THE SINE AND COSINE OF THE ANGLE OF ROTATION |
| 143 |
! |
| 144 |
!Print *, '3 Q=',Q |
| 145 |
!Print *, 'R=',R |
| 146 |
Q = Q - R |
| 147 |
VV = SQRT(4.*P*P+Q*Q) |
| 148 |
!Print *, 'VV=',VV |
| 149 |
!Print *, 'P=',P |
| 150 |
!Print *, 'Q=',Q |
| 151 |
C = SQRT((VV+Q)/(2.*VV)) |
| 152 |
!Print *, 'C=', C |
| 153 |
S = P/(VV*C) |
| 154 |
!Print *, 'S=', S |
| 155 |
! |
| 156 |
! APPLY THE ROTATION TO A |
| 157 |
! |
| 158 |
60 DO 70 I=1,M |
| 159 |
R = A(I,J) |
| 160 |
A(I,J) = R*C + A(I,K)*S |
| 161 |
A(I,K) = -R*S + A(I,K)*C |
| 162 |
70 CONTINUE |
| 163 |
!if (debug) WRITE(*,*) '=======After rotation of A=======' |
| 164 |
! |
| 165 |
! APPLY THE ROTATION TO V |
| 166 |
! |
| 167 |
DO 80 I=1,N |
| 168 |
R = V(I,J) |
| 169 |
V(I,J) = R*C + V(I,K)*S |
| 170 |
!Print *, 'V(I,J)=', V(I,J) |
| 171 |
!Print *, 'R=', R |
| 172 |
V(I,K) = -R*S + V(I,K)*C |
| 173 |
!Print *, 'V(I,K)=', V(I,K) |
| 174 |
80 CONTINUE |
| 175 |
GO TO 100 |
| 176 |
90 ICOUNT = ICOUNT - 1 |
| 177 |
100 CONTINUE |
| 178 |
110 CONTINUE |
| 179 |
!if (debug) WRITE(*,*) '=======After rotation of V=======' |
| 180 |
! |
| 181 |
! PRINT THE NUMBER OF SWEEPS AND ROTATIONS |
| 182 |
! |
| 183 |
IF (NOUT.GT.0) WRITE (NOUT,9997) SWEEP, ICOUNT |
| 184 |
! |
| 185 |
! CHECK NUMBER OF SWEEPS AND ROTATIONS (TERMINATION TEST) |
| 186 |
! |
| 187 |
IF (ICOUNT.GT.0 .AND. SWEEP.LT.LIMIT) GO TO 30 |
| 188 |
IF (SWEEP.GE.LIMIT .AND. NOUT.GT.0) WRITE (NOUT,9996) |
| 189 |
DO 160 J=1,N |
| 190 |
Q = 0. |
| 191 |
DO 120 I=1,M |
| 192 |
Q = Q + A(I,J)*A(I,J) |
| 193 |
120 CONTINUE |
| 194 |
! |
| 195 |
! ARBITRARY RANK DECISION |
| 196 |
! |
| 197 |
IF (J.EQ.1) VV = 1.0E-3*SQRT(Q) |
| 198 |
IF (SQRT(Q).LE.VV) GO TO 140 |
| 199 |
DO 130 I=1,M |
| 200 |
A(I,J) = A(I,J)/Q |
| 201 |
130 CONTINUE |
| 202 |
GO TO 160 |
| 203 |
140 DO 150 I=1,M |
| 204 |
A(I,J) = 0.0 |
| 205 |
150 CONTINUE |
| 206 |
160 CONTINUE |
| 207 |
|
| 208 |
!if (debug) WRITE(*,*) '=======After termination test=======' |
| 209 |
|
| 210 |
DO 190 I=1,N |
| 211 |
DO 180 J=1,M |
| 212 |
TEMP = 0.0 |
| 213 |
DO 170 K=1,N |
| 214 |
TEMP = TEMP + V(I,K)*A(J,K) |
| 215 |
!Print *, 'V(I,K)*A(J,K)=', V(I,K)*A(J,K) |
| 216 |
!Print *, 'V(I,K)=', V(I,K) |
| 217 |
!Print *, 'A(J,K)=', A(J,K) |
| 218 |
170 CONTINUE |
| 219 |
AIN(I,J) = TEMP |
| 220 |
!Print *, 'AIN(I,J)=', AIN(I,J) |
| 221 |
180 CONTINUE |
| 222 |
190 CONTINUE |
| 223 |
RETURN |
| 224 |
|
| 225 |
!Print *, 'AIN=' |
| 226 |
!DO J=1,NAIN |
| 227 |
!WRITE(*,'(50(1X,F10.2))') AIN(:,J) |
| 228 |
!Print *, AIN(:,J) |
| 229 |
!END DO |
| 230 |
!Print *, 'End of AIN' |
| 231 |
|
| 232 |
! |
| 233 |
! IF AN ERROR OCCURED THE APPROPRIATE MESSAGE IS PRINTED |
| 234 |
! AND FAIL IS SET TO .TRUE. |
| 235 |
! |
| 236 |
200 IF (NOUT.GT.0) WRITE (NOUT,9995) |
| 237 |
FAIL = .TRUE. |
| 238 |
RETURN |
| 239 |
210 IF (NOUT.GT.0) WRITE (NOUT,9999) N, NA |
| 240 |
FAIL = .TRUE. |
| 241 |
RETURN |
| 242 |
220 IF (NOUT.GT.0) WRITE (NOUT,9998) M, MA |
| 243 |
FAIL = .TRUE. |
| 244 |
RETURN |
| 245 |
230 IF (NOUT.GT.0) WRITE (NOUT,9994) M, MA, N, NA |
| 246 |
FAIL = .TRUE. |
| 247 |
RETURN |
| 248 |
! |
| 249 |
! * * * * * * * * * * * * * * * * * * * |
| 250 |
! |
| 251 |
9999 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR N IN GINVSE NA= , I5, 5(/)) |
| 252 |
9998 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR M IN GINVSE MA= , I5, 5(/)) |
| 253 |
9997 FORMAT (8H SWEEP =, I4, 2H , I4, 19H JACOBI ROTATIONS P,8HERFORMED) |
| 254 |
9996 FORMAT (21H SWEEP LIMIT REACHED ) |
| 255 |
9995 FORMAT (38H0**ERROR** N GREATER THAN M IS INVALID,21H AND GINVSE WILL FAIL, 5(/)) |
| 256 |
9994 FORMAT (46H0**ERROR** SOME DIMENSIONS ARE LESS THAN TWO I,8HN GINVSE & |
| 257 |
/10X, 3HM =, I5, 5H MA =,I5, 4H N =, I5,5H NA =, I5, 5(/)) |
| 258 |
! |
| 259 |
! * * * END OF SUBROUTINE GINVSE * * * |
| 260 |
! |
| 261 |
END SUBROUTINE GINVSE |
| 262 |
|