]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/02/15 17:48:48 mclareni | |
6 | * Kernlib | |
7 | * | |
8 | * | |
9 | #include "kernnum/pilot.h" | |
10 | SUBROUTINE DEQINV(N,A,IDIM,R,IFAIL,K,B) | |
11 | REAL R(N),T1,T2,T3 | |
12 | DOUBLE PRECISION A(IDIM,N),B(IDIM,K),DET,TEMP,S, | |
13 | $ B1,B2,C11,C12,C13,C21,C22,C23,C31,C32,C33 | |
14 | CHARACTER*6 NAME | |
15 | DATA NAME/'DEQINV'/,KPRNT/1/ | |
16 | C | |
17 | C ****************************************************************** | |
18 | C | |
19 | C REPLACES B BY THE SOLUTION X OF A*X=B, AND REPLACES A BY ITS IN- | |
20 | C VERSE. | |
21 | C | |
22 | C N ORDER OF THE SQUARE MATRIX IN ARRAY A. | |
23 | C | |
24 | C A (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY CONTAINING | |
25 | C AN N BY N MATRIX. | |
26 | C | |
27 | C IDIM FIRST DIMENSION PARAMETER OF ARRAYS A AND B. | |
28 | C | |
29 | C R (REAL) WORKING VECTOR OF LENGTH NOT LESS THAN N. | |
30 | C | |
31 | C IFAIL OUTPUT PARAMETER. IFAIL= 0 ... NORMAL EXIT. | |
32 | C IFAIL=-1 ... SINGULAR MATRIX. | |
33 | C | |
34 | C K NUMBER OF COLUMNS OF THE MATRIX IN ARRAY B. | |
35 | C | |
36 | C B (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY CONTAINING | |
37 | C AN N BY K MATRIX. | |
38 | C | |
39 | C CALLS ... DFACT, DFINV, F010PR, ABEND. | |
40 | C | |
41 | C ****************************************************************** | |
42 | C | |
43 | C TEST FOR PARAMETER ERRORS. | |
44 | C | |
45 | IF((N.LT.1).OR.(N.GT.IDIM).OR.(K.LT.1)) GO TO 10 | |
46 | C | |
47 | C TEST FOR N.LE.3. | |
48 | C | |
49 | IF(N.GT.3) GO TO 9 | |
50 | IFAIL=0 | |
51 | IF(N.LT.3) GO TO 5 | |
52 | C | |
53 | C N=3 CASE. | |
54 | C | |
55 | C COMPUTE COFACTORS. | |
56 | C11=A(2,2)*A(3,3)-A(2,3)*A(3,2) | |
57 | C12=A(2,3)*A(3,1)-A(2,1)*A(3,3) | |
58 | C13=A(2,1)*A(3,2)-A(2,2)*A(3,1) | |
59 | C21=A(3,2)*A(1,3)-A(3,3)*A(1,2) | |
60 | C22=A(3,3)*A(1,1)-A(3,1)*A(1,3) | |
61 | C23=A(3,1)*A(1,2)-A(3,2)*A(1,1) | |
62 | C31=A(1,2)*A(2,3)-A(1,3)*A(2,2) | |
63 | C32=A(1,3)*A(2,1)-A(1,1)*A(2,3) | |
64 | C33=A(1,1)*A(2,2)-A(1,2)*A(2,1) | |
65 | T1=ABS(SNGL(A(1,1))) | |
66 | T2=ABS(SNGL(A(2,1))) | |
67 | T3=ABS(SNGL(A(3,1))) | |
68 | C | |
69 | C (SET TEMP=PIVOT AND DET=PIVOT*DET.) | |
70 | IF(T1.GE.T2) GO TO 1 | |
71 | IF(T3.GE.T2) GO TO 2 | |
72 | C (PIVOT IS A21) | |
73 | TEMP=A(2,1) | |
74 | DET=C13*C32-C12*C33 | |
75 | GO TO 3 | |
76 | 1 IF(T3.GE.T1) GO TO 2 | |
77 | C (PIVOT IS A11) | |
78 | TEMP=A(1,1) | |
79 | DET=C22*C33-C23*C32 | |
80 | GO TO 3 | |
81 | C (PIVOT IS A31) | |
82 | 2 TEMP=A(3,1) | |
83 | DET=C23*C12-C22*C13 | |
84 | C | |
85 | C SET ELEMENTS OF INVERSE IN A. | |
86 | 3 IF(DET.EQ.0D0) GO TO 11 | |
87 | S=TEMP/DET | |
88 | A(1,1)=S*C11 | |
89 | A(1,2)=S*C21 | |
90 | A(1,3)=S*C31 | |
91 | A(2,1)=S*C12 | |
92 | A(2,2)=S*C22 | |
93 | A(2,3)=S*C32 | |
94 | A(3,1)=S*C13 | |
95 | A(3,2)=S*C23 | |
96 | A(3,3)=S*C33 | |
97 | C | |
98 | C REPLACE B BY AINV*B. | |
99 | DO 4 J=1,K | |
100 | B1=B(1,J) | |
101 | B2=B(2,J) | |
102 | B(1,J)=A(1,1)*B1+A(1,2)*B2+A(1,3)*B(3,J) | |
103 | B(2,J)=A(2,1)*B1+A(2,2)*B2+A(2,3)*B(3,J) | |
104 | B(3,J)=A(3,1)*B1+A(3,2)*B2+A(3,3)*B(3,J) | |
105 | 4 CONTINUE | |
106 | RETURN | |
107 | C | |
108 | 5 IF(N.LT.2) GO TO 7 | |
109 | C | |
110 | C N=2 CASE BY CRAMERS RULE. | |
111 | C | |
112 | DET=A(1,1)*A(2,2)-A(1,2)*A(2,1) | |
113 | IF(DET.EQ.0D0) GO TO 11 | |
114 | S=1D0/DET | |
115 | C11 =S*A(2,2) | |
116 | A(1,2)=-S*A(1,2) | |
117 | A(2,1)=-S*A(2,1) | |
118 | A(2,2)=S*A(1,1) | |
119 | A(1,1)=C11 | |
120 | DO 6 J=1,K | |
121 | B1=B(1,J) | |
122 | B(1,J)=C11*B1+A(1,2)*B(2,J) | |
123 | B(2,J)=A(2,1)*B1+A(2,2)*B(2,J) | |
124 | 6 CONTINUE | |
125 | RETURN | |
126 | C | |
127 | C N=1 CASE. | |
128 | C | |
129 | 7 IF(A(1,1).EQ.0D0) GO TO 11 | |
130 | A(1,1)=1D0/A(1,1) | |
131 | DO 8 J=1,K | |
132 | B(1,J)=A(1,1)*B(1,J) | |
133 | 8 CONTINUE | |
134 | RETURN | |
135 | C | |
136 | C N.GT.3 CASES. FACTORIZE MATRIX, INVERT AND SOLVE SYSTEM. | |
137 | C | |
138 | 9 CALL DFACT(N,A,IDIM,R,IFAIL,DET,JFAIL) | |
139 | IF(IFAIL.NE.0) RETURN | |
140 | CALL DFEQN(N,A,IDIM,R,K,B) | |
141 | CALL DFINV(N,A,IDIM,R) | |
142 | RETURN | |
143 | C | |
144 | C ERROR EXITS. | |
145 | C | |
146 | 10 IFAIL=+1 | |
147 | CALL F010PR(NAME,N,IDIM,K,KPRNT) | |
148 | RETURN | |
149 | C | |
150 | 11 IFAIL=-1 | |
151 | RETURN | |
152 | C | |
153 | END |