home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Frostbyte's 1980s DOS Shareware Collection
/
floppyshareware.zip
/
floppyshareware
/
DOOG
/
PCSSP2.ZIP
/
MATINSL.ZIP
/
DSINV.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
4KB
|
123 lines
C
C ..................................................................
C
C SUBROUTINE DSINV
C
C PURPOSE
C INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
C
C USAGE
C CALL DSINV(A,N,EPS,IER)
C
C DESCRIPTION OF PARAMETERS
C A - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN
C SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT
C MATRIX.
C ON RETURN A CONTAINS THE RESULTANT UPPER
C TRIANGULAR MATRIX IN DOUBLE PRECISION.
C N - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
C EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED
C AS RELATIVE TOLERANCE FOR TEST ON LOSS OF
C SIGNIFICANCE.
C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C IER=0 - NO ERROR
C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
C TER N OR BECAUSE SOME RADICAND IS NON-
C POSITIVE (MATRIX A IS NOT POSITIVE
C DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
C FICANCE)
C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI-
C CANCE. THE RADICAND FORMED AT FACTORIZA-
C TION STEP K+1 WAS STILL POSITIVE BUT NO
C LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
C
C REMARKS
C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
C IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
C LAR MATRIX IS STORED COLUMNWISE TOO.
C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
C CALCULATED RADICANDS ARE POSITIVE.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C DMFSD
C
C METHOD
C SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD.
C
C ..................................................................
C
SUBROUTINE DSINV(A,N,EPS,IER)
C
C
DIMENSION A(1)
DOUBLE PRECISION A,DIN,WORK
C
C FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMFSD
C A = TRANSPOSE(T) * T
CALL DMFSD(A,N,EPS,IER)
IF(IER) 9,1,1
C
C INVERT UPPER TRIANGULAR MATRIX T
C PREPARE INVERSION-LOOP
1 IPIV=N*(N+1)/2
IND=IPIV
C
C INITIALIZE INVERSION-LOOP
DO 6 I=1,N
DIN=1.D0/A(IPIV)
A(IPIV)=DIN
MIN=N
KEND=I-1
LANF=N-KEND
IF(KEND) 5,5,2
2 J=IND
C
C INITIALIZE ROW-LOOP
DO 4 K=1,KEND
WORK=0.D0
MIN=MIN-1
LHOR=IPIV
LVER=J
C
C START INNER LOOP
DO 3 L=LANF,MIN
LVER=LVER+1
LHOR=LHOR+L
3 WORK=WORK+A(LVER)*A(LHOR)
C END OF INNER LOOP
C
A(J)=-WORK*DIN
4 J=J-MIN
C END OF ROW-LOOP
C
5 IPIV=IPIV-MIN
6 IND=IND-1
C END OF INVERSION-LOOP
C
C CALCULATE INVERSE(A) BY MEANS OF INVERSE(T)
C INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T))
C INITIALIZE MULTIPLICATION-LOOP
DO 8 I=1,N
IPIV=IPIV+I
J=IPIV
C
C INITIALIZE ROW-LOOP
DO 8 K=I,N
WORK=0.D0
LHOR=J
C
C START INNER LOOP
DO 7 L=K,N
LVER=LHOR+K-I
WORK=WORK+A(LHOR)*A(LVER)
7 LHOR=LHOR+L
C END OF INNER LOOP
C
A(J)=WORK
8 J=J+K
C END OF ROW- AND MULTIPLICATION-LOOP
C
9 RETURN
END