Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaSmallMatInv.c
Go to the documentation of this file.
1
14#include "fasp.h"
15#include "fasp_functs.h"
16
17#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
19/*---------------------------------*/
20/*-- Public Functions --*/
21/*---------------------------------*/
22
33void fasp_smat_inv_nc2 (REAL *a)
34{
35 const REAL a0 = a[0], a1 = a[1];
36 const REAL a2 = a[2], a3 = a[3];
37
38 const REAL det = a0*a3 - a1*a2;
39
40 if ( ABS(det) < SMALLREAL ) {
41 printf("### WARNING: Matrix is nearly singular, det = %e! Ignore.\n", det);
42 printf("##----------------------------------------------\n");
43 printf("## %12.5e %12.5e \n", a0, a1);
44 printf("## %12.5e %12.5e \n", a2, a3);
45 printf("##----------------------------------------------\n");
46
47 a[0] = 1.0; a[1] = 0.0;
48 a[2] = 0.0; a[3] = 1.0;
49 }
50 else {
51 REAL det_inv = 1.0 / det;
52 a[0] = a3 * det_inv; a[1] = -a1 * det_inv;
53 a[2] = -a2 * det_inv; a[3] = a0 * det_inv;
54 }
55}
56
67void fasp_smat_inv_nc3 (REAL *a)
68{
69 const REAL a0 = a[0], a1 = a[1], a2 = a[2];
70 const REAL a3 = a[3], a4 = a[4], a5 = a[5];
71 const REAL a6 = a[6], a7 = a[7], a8 = a[8];
72
73 const REAL M0 = a4*a8-a5*a7, M3 = a2*a7-a1*a8, M6 = a1*a5-a2*a4;
74 const REAL M1 = a5*a6-a3*a8, M4 = a0*a8-a2*a6, M7 = a2*a3-a0*a5;
75 const REAL M2 = a3*a7-a4*a6, M5 = a1*a6-a0*a7, M8 = a0*a4-a1*a3;
76
77 const REAL det = a0*M0+a3*M3+a6*M6;
78
79 if ( ABS(det) < SMALLREAL ) {
80 printf("### WARNING: Matrix is nearly singular, det = %e! Ignore.\n", det);
81 printf("##----------------------------------------------\n");
82 printf("## %12.5e %12.5e %12.5e \n", a0, a1, a2);
83 printf("## %12.5e %12.5e %12.5e \n", a3, a4, a5);
84 printf("## %12.5e %12.5e %12.5e \n", a6, a7, a8);
85 printf("##----------------------------------------------\n");
86
87 a[0] = 1.0; a[1] = 0.0; a[2] = 0.0;
88 a[3] = 0.0; a[4] = 1.0; a[5] = 0.0;
89 a[6] = 0.0; a[7] = 0.0; a[8] = 1.0;
90 }
91 else {
92 REAL det_inv = 1.0/det;
93 a[0] = M0*det_inv; a[1] = M3*det_inv; a[2] = M6*det_inv;
94 a[3] = M1*det_inv; a[4] = M4*det_inv; a[5] = M7*det_inv;
95 a[6] = M2*det_inv; a[7] = M5*det_inv; a[8] = M8*det_inv;
96 }
97}
98
111void fasp_smat_inv_nc4 (REAL *a)
112{
113 const REAL a11 = a[0], a12 = a[1], a13 = a[2], a14 = a[3];
114 const REAL a21 = a[4], a22 = a[5], a23 = a[6], a24 = a[7];
115 const REAL a31 = a[8], a32 = a[9], a33 = a[10], a34 = a[11];
116 const REAL a41 = a[12], a42 = a[13], a43 = a[14], a44 = a[15];
117
118 const REAL M11 = a22*a33*a44 + a23*a34*a42 + a24*a32*a43 - a22*a34*a43 - a23*a32*a44 - a24*a33*a42;
119 const REAL M12 = a12*a34*a43 + a13*a32*a44 + a14*a33*a42 - a12*a33*a44 - a13*a34*a42 - a14*a32*a43;
120 const REAL M13 = a12*a23*a44 + a13*a24*a42 + a14*a22*a43 - a12*a24*a43 - a13*a22*a44 - a14*a23*a42;
121 const REAL M14 = a12*a24*a33 + a13*a22*a34 + a14*a23*a32 - a12*a23*a34 - a13*a24*a32 - a14*a22*a33;
122 const REAL M21 = a21*a34*a43 + a23*a31*a44 + a24*a33*a41 - a21*a33*a44 - a23*a34*a41 - a24*a31*a43;
123 const REAL M22 = a11*a33*a44 + a13*a34*a41 + a14*a31*a43 - a11*a34*a43 - a13*a31*a44 - a14*a33*a41;
124 const REAL M23 = a11*a24*a43 + a13*a21*a44 + a14*a23*a41 - a11*a23*a44 - a13*a24*a41 - a14*a21*a43;
125 const REAL M24 = a11*a23*a34 + a13*a24*a31 + a14*a21*a33 - a11*a24*a33 - a13*a21*a34 - a14*a23*a31;
126 const REAL M31 = a21*a32*a44 + a22*a34*a41 + a24*a31*a42 - a21*a34*a42 - a22*a31*a44 - a24*a32*a41;
127 const REAL M32 = a11*a34*a42 + a12*a31*a44 + a14*a32*a41 - a11*a32*a44 - a12*a34*a41 - a14*a31*a42;
128 const REAL M33 = a11*a22*a44 + a12*a24*a41 + a14*a21*a42 - a11*a24*a42 - a12*a21*a44 - a14*a22*a41;
129 const REAL M34 = a11*a24*a32 + a12*a21*a34 + a14*a22*a31 - a11*a22*a34 - a12*a24*a31 - a14*a21*a32;
130 const REAL M41 = a21*a33*a42 + a22*a31*a43 + a23*a32*a41 - a21*a32*a43 - a22*a33*a41 - a23*a31*a42;
131 const REAL M42 = a11*a32*a43 + a12*a33*a41 + a13*a31*a42 - a11*a33*a42 - a12*a31*a43 - a13*a32*a41;
132 const REAL M43 = a11*a23*a42 + a12*a21*a43 + a13*a22*a41 - a11*a22*a43 - a12*a23*a41 - a13*a21*a42;
133 const REAL M44 = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31;
134
135 const REAL det = a11*M11 + a12*M21 + a13*M31 + a14*M41;
136
137 if ( ABS(det) < SMALLREAL ) {
138 printf("### WARNING: Matrix is nearly singular, det = %e! Ignore.\n", det);
139 printf("##----------------------------------------------\n");
140 printf("## %12.5e %12.5e %12.5e %12.5e\n", a11, a12, a13, a14);
141 printf("## %12.5e %12.5e %12.5e %12.5e\n", a21, a22, a23, a24);
142 printf("## %12.5e %12.5e %12.5e %12.5e\n", a31, a32, a33, a34);
143 printf("## %12.5e %12.5e %12.5e %12.5e\n", a41, a42, a43, a44);
144 printf("##----------------------------------------------\n");
145
146 a[0] = 1.0; a[1] = 0.0; a[2] = 0.0; a[3] = 0.0;
147 a[4] = 0.0; a[5] = 1.0; a[6] = 0.0; a[7] = 0.0;
148 a[8] = 0.0; a[9] = 0.0; a[10] = 1.0; a[11] = 0.0;
149 a[12] = 0.0; a[13] = 0.0; a[14] = 0.0; a[15] = 1.0;
150 }
151 else {
152 REAL det_inv = 1.0 / det;
153 a[0] = M11 * det_inv; a[1] = M12 * det_inv; a[2] = M13 * det_inv; a[3] = M14 * det_inv;
154 a[4] = M21 * det_inv; a[5] = M22 * det_inv; a[6] = M23 * det_inv; a[7] = M24 * det_inv;
155 a[8] = M31 * det_inv; a[9] = M32 * det_inv; a[10] = M33 * det_inv; a[11] = M34 * det_inv;
156 a[12] = M41 * det_inv; a[13] = M42 * det_inv; a[14] = M43 * det_inv; a[15] = M44 * det_inv;
157 }
158}
159
170void fasp_smat_inv_nc5 (REAL *a)
171{
172 const REAL a0=a[0], a1=a[1], a2=a[2], a3=a[3], a4=a[4];
173 const REAL a5=a[5], a6=a[6], a7=a[7], a8=a[8], a9=a[9];
174 const REAL a10=a[10], a11=a[11], a12=a[12], a13=a[13], a14=a[14];
175 const REAL a15=a[15], a16=a[16], a17=a[17], a18=a[18], a19=a[19];
176 const REAL a20=a[20], a21=a[21], a22=a[22], a23=a[23], a24=a[24];
177
178 REAL det0, det1, det2, det3, det4, det;
179
180 det0 = a6 * ( a12 * (a18*a24-a19*a23) + a17 * (a14*a23-a13*a24) + a22 * (a13*a19 - a14*a18) );
181 det0 += a11 * ( a7 * (a19*a23-a18*a24) + a17 * (a8*a24 -a9*a23 ) + a22 * (a9*a18 - a8*a19) );
182 det0 += a16 * ( a7 * (a13*a24-a14*a23) + a12 * (a9*a23 -a8*a24 ) + a22 * (a8*a14 - a9*a13) );
183 det0 += a21 * ( a17 * (a9*a13 -a8*a14 ) + a7 * (a14*a18-a13*a19) + a12 * (a8*a19 - a9*a18) );
184
185 det1 = a1 * ( a22 * (a14*a18-a13*a19) + a12 * (a19*a23-a18*a24) + a17 * (a13*a24 - a14*a23) );
186 det1 += a11 * ( a17 * (a4*a23 - a3*a24) + a2 * (a18*a24-a19*a23) + a22 * (a3*a19 - a4*a18) );
187 det1 += a16 * ( a12 * (a3*a24 - a4*a23) + a2 * (a14*a23-a13*a24) + a22 * (a4*a13 - a3*a14) );
188 det1 += a21 * ( a2 * (a13*a19-a14*a18) + a12 * (a4*a18 -a3*a19 ) + a17 * (a3*a14 - a4*a13) );
189
190 det2 = a1 * ( a7 * (a18*a24-a19*a23) + a17 * (a9*a23-a8*a24) + a22 * (a8*a19 - a9*a18) );
191 det2 += a6 * ( a2 * (a19*a23-a18*a24) + a17 * (a3*a24-a4*a23) + a22 * (a4*a18 - a3*a19) );
192 det2 += a16 * ( a2 * (a8*a24 -a9*a23 ) + a7 * (a4*a23-a3*a24) + a22 * (a3*a9 - a4*a8) );
193 det2 += a21 * ( a7 * (a3*a19 -a4*a18 ) + a2 * (a9*a18-a8*a19) + a17 * (a4*a8 - a3*a9) );
194
195 det3 = a1 * ( a12* (a8*a24 -a9*a23) + a7 * (a14*a23-a13*a24) + a22 * (a9*a13 - a8*a14) );
196 det3 += a6 * ( a2 * (a13*a24-a14*a23) + a12 * (a4*a23 -a3*a24 ) + a22 * (a3*a14 - a4*a13) );
197 det3 += a11 * ( a7 * (a3*a24 -a4*a23) + a2 * (a9*a23 -a8*a24 ) + a22 * (a4*a8 - a3*a9) );
198 det3 += a21 * ( a2 * (a8*a14 -a9*a13) + a7 * (a4*a13 -a3*a14 ) + a12 * (a3*a9 - a4*a8) );
199
200 det4 = a1 * ( a7 * (a13*a19-a14*a18) + a12 * (a9*a18 -a8*a19 ) + a17 * (a8*a14 - a9*a13) );
201 det4 += a6 * ( a12* (a3*a19 -a4*a18 ) + a17 * (a4*a13 -a3*a14 ) + a2 * (a14*a18- a13*a19));
202 det4 += a11 * ( a2 * (a8*a19 -a9*a18 ) + a7 * (a4*a18 -a3*a19 ) + a17 * (a3*a9 - a4*a8) );
203 det4 += a16 * ( a7 * (a3*a14 -a4*a13 ) + a2 * (a9*a13 -a8*a14 ) + a12 * (a4*a8 - a3*a9) );
204
205 det = det0*a0 + det1*a5+ det2*a10 + det3*a15 + det4*a20;
206
207 if ( ABS(det) < SMALLREAL ) {
208 printf("### WARNING: Matrix is nearly singular, det = %e! Ignore.\n", det);
209 printf("##----------------------------------------------\n");
210 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a0, a1, a2, a3, a4);
211 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a5, a6, a7, a8, a9);
212 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a10, a11, a12, a13, a14);
213 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a15, a16, a17, a18, a19);
214 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a20, a21, a22, a23, a24);
215 printf("##----------------------------------------------\n");
216
217 a[0] = 1.0; a[1] = 0.0; a[2] = 0.0; a[3] = 0.0; a[4] = 0.0;
218 a[5] = 0.0; a[6] = 1.0; a[7] = 0.0; a[8] = 0.0; a[9] = 0.0;
219 a[10] = 0.0; a[11] = 0.0; a[12] = 1.0; a[13] = 0.0; a[14] = 0.0;
220 a[15] = 0.0; a[16] = 0.0; a[17] = 0.0; a[18] = 1.0; a[19] = 0.0;
221 a[20] = 0.0; a[21] = 0.0; a[22] = 0.0; a[23] = 0.0; a[24] = 1.0;
222 }
223 else {
224 REAL det_inv = 1 / det;
225
226 a[0] = a6 * (a12 * a18 * a24 - a12 * a19 * a23 - a17 * a13 * a24 + a17 * a14 * a23 + a22 * a13 * a19 - a22 * a14 * a18);
227 a[0] += a11 * (a7 * a19 * a23 - a7 * a18 * a24 + a17 * a8 * a24 - a17 * a9 * a23 - a22 * a8 * a19 + a22 * a9 * a18);
228 a[0] += a16 * (a7 * a13 * a24 - a7 * a14 * a23 - a12 * a8 * a24 + a12 * a9 * a23 + a22 * a8 * a14 - a22 * a9 * a13);
229 a[0] += a21 * (a7 * a14 * a18 - a7 * a13 * a19 + a12 * a8 * a19 - a12 * a9 * a18 - a17 * a8 * a14 + a17 * a9 * a13);
230 a[0] *= det_inv;
231
232 a[1] = a1 * (a12 * a19 * a23 - a12 * a18 * a24 + a22 * a14 * a18 - a17 * a14 * a23 - a22 * a13 * a19 + a17 * a13 * a24);
233 a[1] += a11 * (a22 * a3 * a19 + a2 * a18 * a24 - a17 * a3 * a24 - a22 * a4 * a18 - a2 * a19 * a23 + a17 * a4 * a23);
234 a[1] += a16 * (a12 * a3 * a24 - a12 * a4 * a23 - a22 * a3 * a14 + a2 * a14 * a23 + a22 * a4 * a13 - a2 * a13 * a24);
235 a[1] += a21 * (a12 * a4 * a18 - a12 * a3 * a19 - a2 * a14 * a18 - a17 * a4 * a13 + a2 * a13 * a19 + a17 * a3 * a14);
236 a[1] *= det_inv;
237
238 a[2] = a1 * (a7 * a18 * a24 - a7 * a19 * a23 - a17 * a8 * a24 + a17 * a9 * a23 + a22 * a8 * a19 - a22 * a9 * a18);
239 a[2] += a6 * (a2 * a19 * a23 - a2 * a18 * a24 + a17 * a3 * a24 - a17 * a4 * a23 - a22 * a3 * a19 + a22 * a4 * a18);
240 a[2] += a16 * (a2 * a8 * a24 - a2 * a9 * a23 - a7 * a3 * a24 + a7 * a4 * a23 + a22 * a3 * a9 - a22 * a4 * a8);
241 a[2] += a21 * (a2 * a9 * a18 - a2 * a8 * a19 + a7 * a3 * a19 - a7 * a4 * a18 - a17 * a3 * a9 + a17 * a4 * a8);
242 a[2] *= det_inv;
243
244 a[3] = a1 * (a12 * a8 * a24 - a12 * a9 * a23 + a7 * a14 * a23 - a7 * a13 * a24 + a22 * a9 * a13 - a22 * a8 * a14);
245 a[3] += a6 * (a12 * a4 * a23 - a12 * a3 * a24 + a22 * a3 * a14 - a22 * a4 * a13 + a2 * a13 * a24 - a2 * a14 * a23);
246 a[3] += a11 * (a7 * a3 * a24 - a7 * a4 * a23 + a22 * a4 * a8 - a22 * a3 * a9 + a2 * a9 * a23 - a2 * a8 * a24);
247 a[3] += a21 * (a12 * a3 * a9 - a12 * a4 * a8 + a2 * a8 * a14 - a2 * a9 * a13 + a7 * a4 * a13 - a7 * a3 * a14);
248 a[3] *= det_inv;
249
250 a[4] = a1 * (a7 * a13 * a19 - a7 * a14 * a18 - a12 * a8 * a19 + a12 * a9 * a18 + a17 * a8 * a14 - a17 * a9 * a13);
251 a[4] += a6 * (a2 * a14 * a18 - a2 * a13 * a19 + a12 * a3 * a19 - a12 * a4 * a18 - a17 * a3 * a14 + a17 * a4 * a13);
252 a[4] += a11 * (a2 * a8 * a19 - a2 * a9 * a18 - a7 * a3 * a19 + a7 * a4 * a18 + a17 * a3 * a9 - a17 * a4 * a8);
253 a[4] += a16 * (a2 * a9 * a13 - a2 * a8 * a14 + a7 * a3 * a14 - a7 * a4 * a13 - a12 * a3 * a9 + a12 * a4 * a8);
254 a[4] *= det_inv;
255
256 a[5] = a5 * (a12 * a19 * a23 - a12 * a18 * a24 + a22 * a14 * a18 - a22 * a13 * a19 + a17 * a13 * a24 - a17 * a14 * a23);
257 a[5] += a20 * (a12 * a9 * a18 - a12 * a8 * a19 + a7 * a13 * a19 - a18 * a7 * a14 + a17 * a8 * a14 - a9 * a17 * a13);
258 a[5] += a15 * (a22 * a9 * a13 - a12 * a9 * a23 + a12 * a24 * a8 + a7 * a14 * a23 - a24 * a7 * a13 - a22 * a14 * a8);
259 a[5] += a10 * (a18 * a7 * a24 - a18 * a22 * a9 - a17 * a8 * a24 + a17 * a9 * a23 + a22 * a8 * a19 - a19 * a23 * a7);
260 a[5] *= det_inv;
261
262 a[6] = a2 * (a19 * a23 * a10 - a14 * a23 * a15 - a18 * a24 * a10 + a18 * a14 * a20 - a13 * a19 * a20 + a24 * a13 * a15);
263 a[6] += a12 * (a18 * a0 * a24 - a18 * a20 * a4 + a3 * a19 * a20 - a19 * a23 * a0 + a4 * a23 * a15 - a24 * a15 * a3);
264 a[6] += a17 * (a4 * a13 * a20 - a13 * a24 * a0 + a14 * a23 * a0 - a3 * a14 * a20 + a24 * a3 * a10 - a4 * a23 * a10);
265 a[6] += a22 * (a14 * a15 * a3 - a18 * a14 * a0 + a18 * a4 * a10 - a4 * a13 * a15 + a13 * a19 * a0 - a3 * a19 * a10);
266 a[6] *= det_inv;
267
268 a[7] = a0 * (a18 * a9 * a22 - a18 * a24 * a7 + a19 * a23 * a7 - a9 * a23 * a17 + a24 * a8 * a17 - a8 * a19 * a22);
269 a[7] += a5 * (a2 * a18 * a24 - a2 * a19 * a23 + a17 * a4 * a23 - a17 * a3 * a24 + a22 * a3 * a19 - a22 * a4 * a18);
270 a[7] += a15 * (a4 * a8 * a22 - a3 * a9 * a22 - a24 * a8 * a2 + a9 * a23 * a2 - a4 * a23 * a7 + a24 * a3 * a7);
271 a[7] += a20 * (a18 * a4 * a7 - a18 * a9 * a2 + a9 * a3 * a17 - a4 * a8 * a17 + a8 * a19 * a2 - a3 * a19 * a7);
272 a[7] *= det_inv;
273
274 a[8] = a0 * (a12 * a9 * a23 - a12 * a24 * a8 + a22 * a14 * a8 - a7 * a14 * a23 + a24 * a7 * a13 - a9 * a22 * a13);
275 a[8] += a5 * (a12 * a3 * a24 - a12 * a4 * a23 - a22 * a3 * a14 + a2 * a14 * a23 - a2 * a13 * a24 + a22 * a4 * a13);
276 a[8] += a10 * (a22 * a9 * a3 - a4 * a22 * a8 + a4 * a7 * a23 - a2 * a9 * a23 + a24 * a2 * a8 - a7 * a24 * a3);
277 a[8] += a20 * (a7 * a14 * a3 - a4 * a7 * a13 + a9 * a2 * a13 + a12 * a4 * a8 - a12 * a9 * a3 - a2 * a14 * a8);
278 a[8] *= det_inv;
279
280 a[9] = a0 * (a12 * a8 * a19 - a12 * a18 * a9 + a18 * a7 * a14 - a8 * a17 * a14 + a17 * a13 * a9 - a7 * a13 * a19);
281 a[9] += a5 * (a2 * a13 * a19 - a2 * a14 * a18 - a12 * a3 * a19 + a12 * a4 * a18 + a17 * a3 * a14 - a17 * a4 * a13);
282 a[9] += a10 * (a18 * a2 * a9 - a18 * a7 * a4 + a3 * a7 * a19 - a2 * a8 * a19 + a17 * a8 * a4 - a3 * a17 * a9);
283 a[9] += a15 * (a8 * a2 * a14 - a12 * a8 * a4 + a12 * a3 * a9 - a3 * a7 * a14 + a7 * a13 * a4 - a2 * a13 * a9);
284 a[9] *= det_inv;
285
286 a[10] = a5 * (a18 * a24 * a11 - a24 * a13 * a16 + a14 * a23 * a16 - a19 * a23 * a11 + a13 * a19 * a21 - a18 * a14 * a21);
287 a[10] += a10 * (a19 * a23 * a6 - a9 * a23 * a16 + a24 * a8 * a16 - a8 * a19 * a21 + a18 * a9 * a21 - a18 * a24 * a6);
288 a[10] += a15 * (a24 * a13 * a6 - a14 * a23 * a6 - a24 * a8 * a11 + a9 * a23 * a11 + a14 * a8 * a21 - a13 * a9 * a21);
289 a[10] += a20 * (a18 * a14 * a6 - a18 * a9 * a11 + a8 * a19 * a11 - a13 * a19 * a6 + a9 * a13 * a16 - a14 * a8 * a16);
290 a[10] *= det_inv;
291
292 a[11] = a4 * (a21 * a13 * a15 - a11 * a23 * a15 + a16 * a23 * a10 - a13 * a16 * a20 + a18 * a11 * a20 - a18 * a21 * a10);
293 a[11] += a14 * (a18 * a0 * a21 - a1 * a18 * a20 + a16 * a3 * a20 - a23 * a0 * a16 + a1 * a23 * a15 - a21 * a3 * a15);
294 a[11] += a19 * (a1 * a13 * a20 - a1 * a23 * a10 + a23 * a0 * a11 + a21 * a3 * a10 - a11 * a3 * a20 - a13 * a0 * a21);
295 a[11] += a24 * (a13 * a0 * a16 - a18 * a0 * a11 + a11 * a3 * a15 + a1 * a18 * a10 - a1 * a13 * a15 - a16 * a3 * a10);
296 a[11] *= det_inv;
297
298 a[12] = a4 * (a5 * a21 * a18 - a18 * a20 * a6 + a20 * a16 * a8 - a5 * a16 * a23 + a15 * a6 * a23 - a21 * a15 * a8);
299 a[12] += a9 * (a1 * a20 * a18 - a1 * a15 * a23 + a0 * a16 * a23 - a18 * a0 * a21 - a20 * a16 * a3 + a15 * a21 * a3);
300 a[12] += a19 * (a20 * a6 * a3 - a5 * a21 * a3 + a0 * a21 * a8 - a23 * a0 * a6 + a1 * a5 * a23 - a1 * a20 * a8);
301 a[12] += a24 * (a1 * a15 * a8 - a0 * a16 * a8 + a18 * a0 * a6 - a1 * a5 * a18 + a5 * a16 * a3 - a6 * a15 * a3);
302 a[12] *= det_inv;
303
304 a[13] = a0 * (a24 * a11 * a8 - a6 * a24 * a13 + a21 * a9 * a13 - a11 * a9 * a23 + a14 * a6 * a23 - a14 * a21 * a8);
305 a[13] += a1 * (a5 * a13 * a24 - a5 * a14 * a23 + a14 * a20 * a8 + a10 * a9 * a23 - a24 * a10 * a8 - a20 * a9 * a13);
306 a[13] += a3 * (a6 * a10 * a24 - a10 * a9 * a21 + a5 * a14 * a21 - a5 * a24 * a11 + a20 * a9 * a11 - a14 * a6 * a20);
307 a[13] += a4 * (a5 * a11 * a23 - a5 * a21 * a13 + a21 * a10 * a8 - a6 * a10 * a23 + a20 * a6 * a13 - a11 * a20 * a8);
308 a[13] *= det_inv;
309
310 a[14] = a0 * (a13 * a19 * a6 - a14 * a18 * a6 - a11 * a19 * a8 + a14 * a16 * a8 + a11 * a18 * a9 - a13 * a16 * a9);
311 a[14] += a1 * (a14 * a18 * a5 - a13 * a19 * a5 + a10 * a19 * a8 - a14 * a15 * a8 - a10 * a18 * a9 + a13 * a15 * a9);
312 a[14] += a3 * (a11 * a19 * a5 - a11 * a15 * a9 + a10 * a16 * a9 - a10 * a19 * a6 + a14 * a15 * a6 - a14 * a16 * a5);
313 a[14] += a4 * (a11 * a15 * a8 - a11 * a18 * a5 + a13 * a16 * a5 - a13 * a15 * a6 + a10 * a18 * a6 - a10 * a16 * a8);
314 a[14] *= det_inv;
315
316 a[15] = a5 * (a19 * a22 * a11 - a24 * a17 * a11 + a12 * a24 * a16 - a22 * a14 * a16 - a12 * a19 * a21 + a17 * a14 * a21);
317 a[15] += a10 * (a24 * a17 * a6 - a19 * a22 * a6 - a24 * a7 * a16 + a22 * a9 * a16 + a19 * a7 * a21 - a17 * a9 * a21);
318 a[15] += a15 * (a22 * a14 * a6 - a9 * a22 * a11 + a24 * a7 * a11 - a12 * a24 * a6 - a7 * a14 * a21 + a12 * a9 * a21);
319 a[15] += a20 * (a12 * a19 * a6 - a17 * a14 * a6 - a19 * a7 * a11 + a9 * a17 * a11 + a7 * a14 * a16 - a12 * a9 * a16);
320 a[15] *= det_inv;
321
322 a[16] = a0 * (a11 * a17 * a24 - a11 * a19 * a22 - a12 * a16 * a24 + a12 * a19 * a21 + a14 * a16 * a22 - a14 * a17 * a21);
323 a[16] += a1 * (a10 * a19 * a22 - a10 * a17 * a24 + a12 * a15 * a24 - a12 * a19 * a20 - a14 * a15 * a22 + a14 * a17 * a20);
324 a[16] += a2 * (a10 * a16 * a24 - a10 * a19 * a21 - a11 * a15 * a24 + a11 * a19 * a20 + a14 * a15 * a21 - a14 * a16 * a20);
325 a[16] += a4 * (a10 * a17 * a21 * +a11 * a15 * a22 - a11 * a17 * a20 - a12 * a15 * a21 + a12 * a16 * a20 - a10 * a16 * a22);
326 a[16] *= det_inv;
327
328 a[17] = a0 * (a21 * a9 * a17 - a6 * a24 * a17 + a19 * a6 * a22 - a0 * a16 * a9 * a22 + a24 * a16 * a7 - a19 * a21 * a7);
329 a[17] += a1 * (a5 * a24 * a17 - a5 * a19 * a22 + a19 * a20 * a7 - a20 * a9 * a17 + a15 * a9 * a22 - a24 * a15 * a7);
330 a[17] += a2 * (a5 * a19 * a21 - a19 * a6 * a20 - a5 * a24 * a16 + a24 * a6 * a15 - a15 * a9 * a21 + a20 * a9 * a16);
331 a[17] += a4 * (a16 * a5 * a22 - a6 * a15 * a22 + a20 * a6 * a17 - a5 * a21 * a17 - a6 * a15 * a22 + a21 * a15 * a7 - a16 * a20 * a7);
332 a[17] *= det_inv;
333
334 a[18] = a0 * (a12 * a24 * a6 - a14 * a22 * a6 - a11 * a24 * a7 + a14 * a21 * a7 + a11 * a22 * a9 - a12 * a21 * a9);
335 a[18] += a1 * (a14 * a22 * a5 - a12 * a24 * a5 + a10 * a24 * a7 - a14 * a20 * a7 - a10 * a22 * a9 + a12 * a20 * a9);
336 a[18] += a2 * (a11 * a24 * a5 - a11 * a20 * a9 + a14 * a20 * a6 - a14 * a21 * a5 + a10 * a21 * a9 - a10 * a24 * a6);
337 a[18] += a4 * (a11 * a20 * a7 - a11 * a22 * a5 + a12 * a21 * a5 + a10 * a22 * a6 - a12 * a20 * a6 - a10 * a21 * a7);
338 a[18] *= det_inv;
339
340 a[19] = a0 * (a12 * a16 * a9 - a6 * a12 * a19 + a6 * a17 * a14 - a17 * a11 * a9 + a11 * a7 * a19 - a16 * a7 * a14);
341 a[19] += a1 * (a5 * a12 * a19 - a5 * a17 * a14 - a12 * a15 * a9 + a17 * a10 * a9 + a15 * a7 * a14 - a10 * a7 * a19);
342 a[19] += a2 * (a11 * a15 * a9 - a5 * a11 * a19 + a5 * a16 * a14 - a6 * a15 * a14 + a6 * a10 * a19 - a16 * a10 * a9);
343 a[19] += a4 * (a5 * a17 * a11 - a5 * a12 * a16 + a12 * a6 * a15 + a10 * a7 * a16 - a17 * a6 * a10 - a15 * a7 * a11);
344 a[19] *= det_inv;
345
346 a[20] = a5 * (a12 * a18 * a21 - a12 * a23 * a16 + a22 * a13 * a16 - a18 * a22 * a11 + a23 * a17 * a11 - a17 * a13 * a21);
347 a[20] += a15 * (a12 * a23 * a6 - a12 * a8 * a21 + a8 * a22 * a11 - a23 * a7 * a11 + a7 * a13 * a21 - a22 * a13 * a6);
348 a[20] += a20 * (a12 * a8 * a16 - a12 * a18 * a6 + a18 * a7 * a11 - a8 * a17 * a11 + a17 * a13 * a6 - a7 * a13 * a16);
349 a[20] += a10 * (a17 * a8 * a21 - a22 * a8 * a16 - a18 * a7 * a21 + a18 * a22 * a6 + a23 * a7 * a16 - a23 * a17 * a6);
350 a[20] *= det_inv;
351
352 a[21] = a0 * (a12 * a23 * a16 - a12 * a18 * a21 + a17 * a13 * a21 + a18 * a22 * a11 - a23 * a17 * a11 - a22 * a13 * a16);
353 a[21] += a1 * (a12 * a18 * a20 - a12 * a23 * a15 + a22 * a13 * a15 + a23 * a17 * a10 - a17 * a13 * a20 - a18 * a22 * a10);
354 a[21] += a2 * (a18 * a21 * a10 - a18 * a11 * a20 - a21 * a13 * a15 + a16 * a13 * a20 - a23 * a16 * a10 + a23 * a11 * a15);
355 a[21] += a3 * (a17 * a11 * a20 - a12 * a16 * a20 + a12 * a21 * a15 - a21 * a17 * a10 - a22 * a11 * a15 + a16 * a22 * a10);
356 a[21] *= det_inv;
357
358 a[22] = a0 * (a18 * a21 * a7 - a18 * a6 * a22 + a23 * a6 * a17 + a16 * a8 * a22 - a21 * a8 * a17 - a23 * a16 * a7);
359 a[22] += a1 * (a5 * a18 * a22 - a5 * a23 * a17 - a15 * a8 * a22 + a20 * a8 * a17 - a18 * a20 * a7 + a23 * a15 * a7);
360 a[22] += a3 * (a16 * a20 * a7 + a6 * a15 * a22 - a6 * a20 * a17 - a5 * a16 * a22 + a5 * a21 * a17 - a21 * a15 * a7);
361 a[22] += a2 * (a5 * a23 * a16 - a5 * a18 * a21 + a18 * a6 * a20 + a15 * a8 * a21 - a20 * a8 * a16 - a23 * a6 * a15);
362 a[22] *= det_inv;
363
364 a[23] = a0 * (a12 * a21 * a8 - a22 * a11 * a8 + a11 * a7 * a23 - a6 * a12 * a23 - a21 * a7 * a13 + a6 * a22 * a13);
365 a[23] += a1 * (a5 * a12 * a23 - a5 * a22 * a13 - a10 * a7 * a23 + a20 * a7 * a13 + a22 * a10 * a8 - a12 * a20 * a8);
366 a[23] += a2 * (a5 * a21 * a13 + a11 * a20 * a8 + a6 * a10 * a23 - a5 * a11 * a23 - a21 * a10 * a8 - a6 * a20 * a13);
367 a[23] += a3 * (a5 * a22 * a11 - a5 * a12 * a21 + a10 * a7 * a21 - a22 * a6 * a10 - a20 * a7 * a11 + a12 * a6 * a20);
368 a[23] *= det_inv;
369
370 a[24] = a0 * (a17 * a11 * a8 - a11 * a7 * a18 + a6 * a12 * a18 - a12 * a16 * a8 + a16 * a7 * a13 - a6 * a17 * a13);
371 a[24] += a1 * (a5 * a17 * a13 - a5 * a12 * a18 + a10 * a7 * a18 + a12 * a15 * a8 - a17 * a10 * a8 - a15 * a7 * a13);
372 a[24] += a2 * (a5 * a11 * a18 - a5 * a16 * a13 + a16 * a10 * a8 + a6 * a15 * a13 - a11 * a15 * a8 - a6 * a10 * a18);
373 a[24] += a3 * (a5 * a12 * a16 + a17 * a6 * a10 - a5 * a17 * a11 - a12 * a6 * a15 - a10 * a7 * a16 + a15 * a7 * a11);
374 a[24] *= det_inv;
375
376 }
377
378 printf("### DEBUG: Check inverse matrix...\n");
379 printf("##----------------------------------------------\n");
380 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n",
381 a0 * a[0] + a1 * a[5] + a2 * a[10] + a3 * a[15] + a4 * a[20],
382 a0 * a[1] + a1 * a[6] + a2 * a[11] + a3 * a[16] + a4 * a[21],
383 a0 * a[2] + a1 * a[7] + a2 * a[12] + a3 * a[17] + a4 * a[22],
384 a0 * a[3] + a1 * a[8] + a2 * a[13] + a3 * a[18] + a4 * a[23],
385 a0 * a[4] + a1 * a[9] + a2 * a[14] + a3 * a[19] + a4 * a[24]);
386 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n",
387 a5 * a[0] + a6 * a[5] + a7 * a[10] + a8 * a[15] + a9 * a[20],
388 a5 * a[1] + a6 * a[6] + a7 * a[11] + a8 * a[16] + a9 * a[21],
389 a5 * a[2] + a6 * a[7] + a7 * a[12] + a8 * a[17] + a9 * a[22],
390 a5 * a[3] + a6 * a[8] + a7 * a[13] + a8 * a[18] + a9 * a[23],
391 a5 * a[4] + a6 * a[9] + a7 * a[14] + a8 * a[19] + a9 * a[24]);
392 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n",
393 a10 * a[0] + a11 * a[5] + a12 * a[10] + a13 * a[15] + a14 * a[20],
394 a10 * a[1] + a11 * a[6] + a12 * a[11] + a13 * a[16] + a14 * a[21],
395 a10 * a[2] + a11 * a[7] + a12 * a[12] + a13 * a[17] + a14 * a[22],
396 a10 * a[3] + a11 * a[8] + a12 * a[13] + a13 * a[18] + a14 * a[23],
397 a10 * a[4] + a11 * a[9] + a12 * a[14] + a13 * a[19] + a14 * a[24]);
398 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n",
399 a15 * a[0] + a16 * a[5] + a17 * a[10] + a18 * a[15] + a19 * a[20],
400 a15 * a[1] + a16 * a[6] + a17 * a[11] + a18 * a[16] + a19 * a[21],
401 a15 * a[2] + a16 * a[7] + a17 * a[12] + a18 * a[17] + a19 * a[22],
402 a15 * a[3] + a16 * a[8] + a17 * a[13] + a18 * a[18] + a19 * a[23],
403 a15 * a[4] + a16 * a[9] + a17 * a[14] + a18 * a[19] + a19 * a[24]);
404 printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n",
405 a20 * a[0] + a21 * a[5] + a22 * a[10] + a23 * a[15] + a24 * a[20],
406 a20 * a[1] + a21 * a[6] + a22 * a[11] + a23 * a[16] + a24 * a[21],
407 a20 * a[2] + a21 * a[7] + a22 * a[12] + a23 * a[17] + a24 * a[22],
408 a20 * a[3] + a21 * a[8] + a22 * a[13] + a23 * a[18] + a24 * a[23],
409 a20 * a[4] + a21 * a[9] + a22 * a[14] + a23 * a[19] + a24 * a[24]);
410 printf("##----------------------------------------------\n");
411}
412
425void fasp_smat_inv_nc7 (REAL *a)
426{
427 fasp_smat_invp_nc(a,7);
428}
429
441void fasp_smat_inv_nc (REAL *a,
442 const INT n)
443{
444 INT i,j,k,l,u,kn,in;
445 REAL alinv;
446
447 for (k=0; k<n; ++k) {
448
449 kn = k*n;
450 l = kn+k;
451
452 if (ABS(a[l]) < SMALLREAL) {
453 printf("### ERROR: Diagonal entry is close to zero! ");
454 printf("diag_%d = %.2e! [%s]\n", k, a[l], __FUNCTION__);
455 exit(ERROR_SOLVER_EXIT);
456 }
457 alinv = 1.0/a[l];
458 a[l] = alinv;
459
460 for (j=0; j<k; ++j) {
461 u = kn+j; a[u] *= alinv;
462 }
463
464 for (j=k+1; j<n; ++j) {
465 u = kn+j; a[u] *= alinv;
466 }
467
468 for (i=0; i<k; ++i) {
469 in = i*n;
470 for (j=0; j<n; ++j)
471 if (j!=k) {
472 u = in+j; a[u] -= a[in+k]*a[kn+j];
473 } // end if (j!=k)
474 }
475
476 for (i=k+1; i<n; ++i) {
477 in = i*n;
478 for (j=0; j<n; ++j)
479 if (j!=k) {
480 u = in+j; a[u] -= a[in+k]*a[kn+j];
481 } // end if (j!=k)
482 }
483
484 for (i=0; i<k; ++i) {
485 u=i*n+k; a[u] *= -alinv;
486 }
487
488 for (i=k+1; i<n; ++i) {
489 u=i*n+k; a[u] *= -alinv;
490 }
491
492 } // end for (k=0; k<n; ++k)
493}
494
508SHORT fasp_smat_invp_nc (REAL *a,
509 const INT n)
510{
511 INT i, j, k, l, ll, u;
512 INT icol = 0, irow = 0;
513 REAL vmax, dum, pivinv, temp;
514
515 INT *work = (INT *)fasp_mem_calloc(3*n,sizeof(INT));
516 INT *indxc = work, *indxr = work+n, *ipiv = work+2*n;
517
518 // ipiv, indxr, and indxc are used for book-keeping on the pivoting.
519 for ( j=0; j<n; j++ ) ipiv[j] = 0;
520
521#if DEBUG_MODE > 1
522 printf("### DEBUG: Matrix block\n");
523 for ( i = 0; i < n; ++i ) {
524 for ( j = 0; j < n; ++j ) {
525 printf(" %10.5e,", a[i * n + j]);
526 }
527 printf("\n");
528 }
529#endif
530
531 // This is the main loop over the columns to be reduced.
532 for ( i=0; i<n; i++ ) {
533
534 // This is the outer loop of the search for a pivot element.
535 vmax = 0.0;
536 for ( j=0; j<n; j++ ) {
537 if ( ipiv[j] != 1 ) {
538 for ( k=0; k<n; k++ ) {
539 if ( ipiv[k] == 0 ) {
540 u = j*n+k;
541 if ( ABS(a[u]) >= vmax ) {
542 vmax = ABS(a[u]); irow = j; icol = k;
543 }
544 }
545 } // end for k
546 }
547 } // end for j
548
549 ++(ipiv[icol]);
550
551 // We now have the pivot element, so we interchange rows, if needed, to put
552 // the pivot element on the diagonal. The columns are not physically
553 // interchanged, only relabeled: indxc[i], the column of the ith pivot
554 // element, is the ith column that is reduced, while indxr[i] is the row in
555 // which that pivot element was originally located. If indxr[i] != indxc[i]
556 // there is an implied column interchange. With this form of bookkeeping,
557 // the inverse matrix will be scrambled by columns.
558 if ( irow != icol ) {
559 for ( l=0; l<n; l++ ) SWAP(a[irow*n+l],a[icol*n+l]);
560 }
561
562 indxr[i] = irow; indxc[i] = icol;
563 u = icol*n+icol;
564 if ( ABS(a[u]) < SMALLREAL ) {
565 printf("### WARNING: The matrix is nearly singular!\n");
566 return ERROR_SOLVER_EXIT;
567 }
568 pivinv = 1.0/a[u]; a[u]=1.0;
569 for ( l=0; l<n; l++ ) a[icol*n+l] *= pivinv;
570
571 for ( ll=0; ll<n; ll++ ) {
572 if ( ll != icol ) {
573 u = ll*n+icol;
574 dum = a[u]; a[u] = 0.0;
575 for ( l=0; l<n; l++ ) a[ll*n+l] -= a[icol*n+l]*dum;
576 }
577 }
578 }
579 // This is the end of the main loop over columns of the reduction.
580
581 // It only remains to unscramble the matrix in view of the column interchanges.
582 for ( l=n-1; l>=0; l-- ) {
583 if ( indxr[l] != indxc[l] )
584 for ( k=0; k<n; k++ ) SWAP(a[k*n+indxr[l]],a[k*n+indxc[l]]);
585 } // And we are done.
586
587 fasp_mem_free(work); work = NULL;
588
589 return FASP_SUCCESS;
590}
591
603SHORT fasp_smat_inv (REAL *a,
604 const INT n)
605{
606 SHORT status = FASP_SUCCESS;
607
608 switch (n) {
609
610 case 2:
611 fasp_smat_inv_nc2(a);
612 break;
613
614 case 3:
615 fasp_smat_inv_nc3(a);
616 break;
617
618 case 4:
619 fasp_smat_inv_nc4(a);
620 break;
621
622 case -5:
623 fasp_smat_inv_nc5(a);
624 break;
625
626 default:
627 status = fasp_smat_invp_nc(a,n);
628 break;
629
630 }
631
632 return status;
633}
634
646REAL fasp_smat_Linf (const REAL *A,
647 const INT n)
648{
649
650 REAL norm = 0.0, value;
651
652 INT i,j;
653
654 for ( i = 0; i < n; i++ ) {
655 for ( value = 0.0, j = 0; j < n; j++ ) {
656 value = value + ABS(A[i*n+j]);
657 }
658 norm = MAX(norm, value);
659 }
660
661 return norm;
662}
663
674void fasp_smat_identity_nc2 (REAL *a)
675{
676 memset(a, 0X0, 4*sizeof(REAL));
677
678 a[0] = 1.0; a[3] = 1.0;
679}
680
691void fasp_smat_identity_nc3 (REAL *a)
692{
693 memset(a, 0X0, 9*sizeof(REAL));
694
695 a[0] = 1.0; a[4] = 1.0; a[8] = 1.0;
696}
697
708void fasp_smat_identity_nc5 (REAL *a)
709{
710 memset(a, 0X0, 25*sizeof(REAL));
711
712 a[0] = 1.0;
713 a[6] = 1.0;
714 a[12] = 1.0;
715 a[18] = 1.0;
716 a[24] = 1.0;
717}
718
729void fasp_smat_identity_nc7 (REAL *a)
730{
731 memset(a, 0X0, 49*sizeof(REAL));
732
733 a[0] = 1.0;
734 a[8] = 1.0;
735 a[16] = 1.0;
736 a[24] = 1.0;
737 a[32] = 1.0;
738 a[40] = 1.0;
739 a[48] = 1.0;
740}
741
754void fasp_smat_identity (REAL *a,
755 const INT n,
756 const INT n2)
757{
758 memset(a, 0X0, n2*sizeof(REAL));
759
760 switch (n) {
761
762 case 2: {
763 a[0] = 1.0;
764 a[3] = 1.0;
765 }
766 break;
767
768 case 3: {
769 a[0] = 1.0;
770 a[4] = 1.0;
771 a[8] = 1.0;
772 }
773 break;
774
775 case 4: {
776 a[0] = 1.0;
777 a[5] = 1.0;
778 a[10] = 1.0;
779 a[15] = 1.0;
780 }
781 break;
782
783 case 5: {
784 a[0] = 1.0;
785 a[6] = 1.0;
786 a[12] = 1.0;
787 a[18] = 1.0;
788 a[24] = 1.0;
789 }
790 break;
791
792 case 6: {
793 a[0] = 1.0;
794 a[7] = 1.0;
795 a[14] = 1.0;
796 a[21] = 1.0;
797 a[28] = 1.0;
798 a[35] = 1.0;
799 }
800 break;
801
802 case 7: {
803 a[0] = 1.0;
804 a[8] = 1.0;
805 a[16] = 1.0;
806 a[24] = 1.0;
807 a[32] = 1.0;
808 a[40] = 1.0;
809 a[48] = 1.0;
810 }
811 break;
812
813 default: {
814 INT l;
815 for (l = 0; l < n; l ++) a[l*n+l] = 1.0;
816 }
817 break;
818 }
819
820}
821
822/*---------------------------------*/
823/*-- End of File --*/
824/*---------------------------------*/
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_smat_Linf
REAL fasp_smat_Linf(const REAL *A, const INT n)
Compute the L infinity norm of A.
Definition: BlaSmallMatInv.c:646
fasp_smat_inv_nc
void fasp_smat_inv_nc(REAL *a, const INT n)
Compute the inverse of a matrix using Gauss Elimination.
Definition: BlaSmallMatInv.c:441
fasp_smat_inv_nc3
void fasp_smat_inv_nc3(REAL *a)
Compute the inverse matrix of a 3*3 full matrix A (in place)
Definition: BlaSmallMatInv.c:67
fasp_smat_identity_nc7
void fasp_smat_identity_nc7(REAL *a)
Set a 7*7 full matrix to be a identity.
Definition: BlaSmallMatInv.c:729
fasp_smat_inv_nc4
void fasp_smat_inv_nc4(REAL *a)
Compute the inverse matrix of a 4*4 full matrix A (in place)
Definition: BlaSmallMatInv.c:111
fasp_smat_identity_nc3
void fasp_smat_identity_nc3(REAL *a)
Set a 3*3 full matrix to be a identity.
Definition: BlaSmallMatInv.c:691
fasp_smat_invp_nc
SHORT fasp_smat_invp_nc(REAL *a, const INT n)
Compute the inverse of a matrix using Gauss Elimination with Pivoting.
Definition: BlaSmallMatInv.c:508
SWAP
#define SWAP(a, b)
Definition: BlaSmallMatInv.c:17
fasp_smat_inv_nc7
void fasp_smat_inv_nc7(REAL *a)
Compute the inverse matrix of a 7*7 matrix a.
Definition: BlaSmallMatInv.c:425
fasp_smat_inv_nc2
void fasp_smat_inv_nc2(REAL *a)
Compute the inverse matrix of a 2*2 full matrix A (in place)
Definition: BlaSmallMatInv.c:33
fasp_smat_inv
SHORT fasp_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: BlaSmallMatInv.c:603
fasp_smat_identity_nc2
void fasp_smat_identity_nc2(REAL *a)
Set a 2*2 full matrix to be a identity.
Definition: BlaSmallMatInv.c:674
fasp_smat_inv_nc5
void fasp_smat_inv_nc5(REAL *a)
Compute the inverse matrix of a 5*5 full matrix A (in place)
Definition: BlaSmallMatInv.c:170
fasp_smat_identity_nc5
void fasp_smat_identity_nc5(REAL *a)
Set a 5*5 full matrix to be a identity.
Definition: BlaSmallMatInv.c:708
fasp_smat_identity
void fasp_smat_identity(REAL *a, const INT n, const INT n2)
Set a n*n full matrix to be a identity.
Definition: BlaSmallMatInv.c:754
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
ERROR_SOLVER_EXIT
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:49
SMALLREAL
#define SMALLREAL
Definition: fasp_const.h:256