Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
smat.c
Go to the documentation of this file.
1 
6 #include "fasp.h"
7 #include "fasp_functs.h"
8 
9 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
11 /*---------------------------------*/
12 /*-- Public Functions --*/
13 /*---------------------------------*/
14 
26 {
27  const REAL a0 = a[0], a1 = a[1];
28  const REAL a2 = a[2], a3 = a[3];
29 
30  const REAL det = a0*a3 - a1*a2;
31  REAL det_inv;
32 
33  if ( ABS(det) < SMALLREAL ) {
34  printf("### WARNING: Matrix is nearly singular! det = %e\n", det);
35  /*
36  printf("##----------------------------------------------\n");
37  printf("## %12.5e %12.5e \n", a0, a1);
38  printf("## %12.5e %12.5e \n", a2, a3);
39  printf("##----------------------------------------------\n");
40  getchar();
41  */
42  }
43 
44  det_inv = 1.0/det;
45 
46  a[0] = a3*det_inv; a[1] = - a1*det_inv;
47 
48  a[2] = - a2*det_inv; a[3] = a0*det_inv;
49 }
50 
62 {
63  const REAL a0 = a[0], a1 = a[1], a2 = a[2];
64  const REAL a3 = a[3], a4 = a[4], a5 = a[5];
65  const REAL a6 = a[6], a7 = a[7], a8 = a[8];
66 
67  const REAL M0 = a4*a8-a5*a7, M3 = a2*a7-a1*a8, M6 = a1*a5-a2*a4;
68  const REAL M1 = a5*a6-a3*a8, M4 = a0*a8-a2*a6, M7 = a2*a3-a0*a5;
69  const REAL M2 = a3*a7-a4*a6, M5 = a1*a6-a0*a7, M8 = a0*a4-a1*a3;
70  const REAL det = a0*M0+a3*M3+a6*M6;
71  REAL det_inv;
72 
73  if ( ABS(det) < SMALLREAL ) {
74  printf("### WARNING: Matrix is nearly singular! det = %e\n", det);
75  /*
76  printf("##----------------------------------------------\n");
77  printf("## %12.5e %12.5e %12.5e \n", a0, a1, a2);
78  printf("## %12.5e %12.5e %12.5e \n", a3, a4, a5);
79  printf("## %12.5e %12.5e %12.5e \n", a5, a6, a7);
80  printf("##----------------------------------------------\n");
81  getchar();
82  */
83 
84  a[0] = 1.0; a[1] = 0.0; a[2] = 0.0;
85 
86  a[3] = 0.0; a[4] = 1.0; a[5] = 0.0;
87 
88  a[6] = 0.0; a[7] = 0.0; a[8] = 1.0;
89 
90  }
91  else {
92 
93  det_inv = 1.0/det;
94 
95  a[0] = M0*det_inv; a[1] = M3*det_inv; a[2] = M6*det_inv;
96 
97  a[3] = M1*det_inv; a[4] = M4*det_inv; a[5] = M7*det_inv;
98 
99  a[6] = M2*det_inv; a[7] = M5*det_inv; a[8] = M8*det_inv;
100  }
101 }
102 
116 {
117  const REAL a11 = a[0], a12 = a[1], a13 = a[2], a14 = a[3];
118  const REAL a21 = a[4], a22 = a[5], a23 = a[6], a24 = a[7];
119  const REAL a31 = a[8], a32 = a[9], a33 = a[10], a34 = a[11];
120  const REAL a41 = a[12], a42 = a[13], a43 = a[14], a44 = a[15];
121 
122  const REAL M11 = a22*a33*a44 + a23*a34*a42 + a24*a32*a43 - a22*a34*a43 - a23*a32*a44 - a24*a33*a42;
123  const REAL M12 = a12*a34*a43 + a13*a32*a44 + a14*a33*a42 - a12*a33*a44 - a13*a34*a42 - a14*a32*a43;
124  const REAL M13 = a12*a23*a44 + a13*a24*a42 + a14*a22*a43 - a12*a24*a43 - a13*a22*a44 - a14*a23*a42;
125  const REAL M14 = a12*a24*a33 + a13*a22*a34 + a14*a23*a32 - a12*a23*a34 - a13*a24*a32 - a14*a22*a33;
126  const REAL M21 = a21*a34*a43 + a23*a31*a44 + a24*a33*a41 - a21*a33*a44 - a23*a34*a41 - a24*a31*a43;
127  const REAL M22 = a11*a33*a44 + a13*a34*a41 + a14*a31*a43 - a11*a34*a43 - a13*a31*a44 - a14*a33*a41;
128  const REAL M23 = a11*a24*a43 + a13*a21*a44 + a14*a23*a41 - a11*a23*a44 - a13*a24*a41 - a14*a21*a43;
129  const REAL M24 = a11*a23*a34 + a13*a24*a31 + a14*a21*a33 - a11*a24*a33 - a13*a21*a34 - a14*a23*a31;
130  const REAL M31 = a21*a32*a44 + a22*a34*a41 + a24*a31*a42 - a21*a34*a42 - a22*a31*a44 - a24*a32*a41;
131  const REAL M32 = a11*a34*a42 + a12*a31*a44 + a14*a32*a41 - a11*a32*a44 - a12*a34*a41 - a14*a31*a42;
132  const REAL M33 = a11*a22*a44 + a12*a24*a41 + a14*a21*a42 - a11*a24*a42 - a12*a21*a44 - a14*a22*a41;
133  const REAL M34 = a11*a24*a32 + a12*a21*a34 + a14*a22*a31 - a11*a22*a34 - a12*a24*a31 - a14*a21*a32;
134  const REAL M41 = a21*a33*a42 + a22*a31*a43 + a23*a32*a41 - a21*a32*a43 - a22*a33*a41 - a23*a31*a42;
135  const REAL M42 = a11*a32*a43 + a12*a33*a41 + a13*a31*a42 - a11*a33*a42 - a12*a31*a43 - a13*a32*a41;
136  const REAL M43 = a11*a23*a42 + a12*a21*a43 + a13*a22*a41 - a11*a22*a43 - a12*a23*a41 - a13*a21*a42;
137  const REAL M44 = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31;
138 
139  const REAL det = a11*M11 + a12*M21 + a13*M31 + a14*M41;
140  REAL det_inv;
141 
142  if ( ABS(det) < SMALLREAL ) {
143  printf("### WARNING: Matrix is nearly singular! det = %e\n", det);
144  /*
145  printf("##----------------------------------------------\n");
146  printf("## %12.5e %12.5e %12.5e \n", a0, a1, a2);
147  printf("## %12.5e %12.5e %12.5e \n", a3, a4, a5);
148  printf("## %12.5e %12.5e %12.5e \n", a5, a6, a7);
149  printf("##----------------------------------------------\n");
150  getchar();
151  */
152  }
153 
154  det_inv = 1.0/det;
155 
156  a[0] = M11*det_inv; a[1] = M12*det_inv; a[2] = M13*det_inv; a[3] = M14*det_inv;
157  a[4] = M21*det_inv; a[5] = M22*det_inv; a[6] = M23*det_inv; a[7] = M24*det_inv;
158  a[8] = M31*det_inv; a[9] = M32*det_inv; a[10] = M33*det_inv; a[11] = M34*det_inv;
159  a[12] = M41*det_inv; a[13] = M42*det_inv; a[14] = M43*det_inv; a[15] = M44*det_inv;
160 
161 }
162 
174 {
175  const REAL a0=a[0], a1=a[1], a2=a[2], a3=a[3], a4=a[4];
176  const REAL a5=a[5], a6=a[6], a7=a[7], a8=a[8], a9=a[9];
177  const REAL a10=a[10], a11=a[11], a12=a[12], a13=a[13], a14=a[14];
178  const REAL a15=a[15], a16=a[16], a17=a[17], a18=a[18], a19=a[19];
179  const REAL a20=a[20], a21=a[21], a22=a[22], a23=a[23], a24=a[24];
180 
181  REAL det0, det1, det2, det3, det4, det, det_inv;
182 
183  det0 = a6 * ( a12 * (a18*a24-a19*a23) + a17 * (a14*a23-a13*a24) + a22 * (a13*a19 - a14*a18) );
184  det0 += a11 * ( a7 * (a19*a23-a18*a24) + a17 * (a8*a24 -a9*a23 ) + a22 * (a9*a18 - a8*a19) );
185  det0 += a16 * ( a7 * (a13*a24-a14*a23) + a12 * (a9*a23 -a8*a24 ) + a22 * (a8*a14 - a9*a13) );
186  det0 += a21 * ( a17 * (a9*a13 -a8*a14 ) + a7 * (a14*a18-a13*a19) + a12 * (a8*a19 - a9*a18) );
187 
188  det1 = a1 * ( a22 * (a14*a18-a13*a19) + a12 * (a19*a23-a18*a24) + a17 * (a13*a24 - a14*a23) );
189  det1 += a11 * ( a17 * (a4*a23 - a3*a24) + a2 * (a18*a24-a19*a23) + a22 * (a3*a19 - a4*a18) );
190  det1 += a16 * ( a12 * (a3*a24 - a4*a23) + a2 * (a14*a23-a13*a24) + a22 * (a4*a13 - a3*a14) );
191  det1 += a21 * ( a2 * (a13*a19-a14*a18) + a12 * (a4*a18 -a3*a19 ) + a17 * (a3*a14 - a4*a13) );
192 
193  det2 = a1 * ( a7 * (a18*a24-a19*a23) + a17 * (a9*a23-a8*a24) + a22 * (a8*a19 - a9*a18) );
194  det2 += a6 * ( a2 * (a19*a23-a18*a24) + a17 * (a3*a24-a4*a23) + a22 * (a4*a18 - a3*a19) );
195  det2 += a16 * ( a2 * (a8*a24 -a9*a23 ) + a7 * (a4*a23-a3*a24) + a22 * (a3*a9 - a4*a8) );
196  det2 += a21 * ( a7 * (a3*a19 -a4*a18 ) + a2 * (a9*a18-a8*a19) + a17 * (a4*a8 - a3*a9) );
197 
198  det3 = a1 * ( a12* (a8*a24 -a9*a23) + a7 * (a14*a23-a13*a24) + a22 * (a9*a13 - a8*a14) );
199  det3 += a6 * ( a2 * (a13*a24-a14*a23) + a12 * (a4*a23 -a3*a24 ) + a22 * (a3*a14 - a4*a13) );
200  det3 += a11 * ( a7 * (a3*a24 -a4*a23) + a2 * (a9*a23 -a8*a24 ) + a22 * (a4*a8 - a3*a9) );
201  det3 += a21 * ( a2 * (a8*a14 -a9*a13) + a7 * (a4*a13 -a3*a14 ) + a12 * (a3*a9 - a4*a8) );
202 
203  det4 = a1 * ( a7 * (a13*a19-a14*a18) + a12 * (a9*a18 -a8*a19 ) + a17 * (a8*a14 - a9*a13) );
204  det4 += a6 * ( a12* (a3*a19 -a4*a18 ) + a17 * (a4*a13 -a3*a14 ) + a2 * (a14*a18- a13*a19));
205  det4 += a11 * ( a2 * (a8*a19 -a9*a18 ) + a7 * (a4*a18 -a3*a19 ) + a17 * (a3*a9 - a4*a8) );
206  det4 += a16 * ( a7 * (a3*a14 -a4*a13 ) + a2 * (a9*a13 -a8*a14 ) + a12 * (a4*a8 - a3*a9) );
207 
208  det = det0*a0 + det1*a5+ det2*a10 + det3*a15 + det4*a20;
209 
210  if ( ABS(det) < SMALLREAL ) {
211  printf("### WARNING: Matrix is nearly singular! det = %e\n", det);
212  /*
213  printf("##----------------------------------------------\n");
214  printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a0, a1, a2, a3, a4);
215  printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a5, a6, a7, a8, a9);
216  printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a10, a11, a12, a13, a14);
217  printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a15, a16, a17, a18, a19);
218  printf("## %12.5e %12.5e %12.5e %12.5e %12.5e\n", a20, a21, a22, a23, a24);
219  printf("##----------------------------------------------\n");
220  getchar();
221  */
222  }
223 
224  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 
390 {
392 }
393 
406  const INT n)
407 {
408  INT i,j,k,l,u,kn,in;
409  REAL alinv;
410 
411  for (k=0; k<n; ++k) {
412 
413  kn = k*n;
414  l = kn+k;
415 
416  if (ABS(a[l]) < SMALLREAL) {
417  printf("### ERROR: Diagonal entry is close to zero! ");
418  printf("diag_%d = %.2e\n", k, a[l]);
419  exit(ERROR_SOLVER_EXIT);
420  }
421  alinv = 1.0/a[l];
422  a[l] = alinv;
423 
424  for (j=0; j<k; ++j) {
425  u = kn+j; a[u] *= alinv;
426  }
427 
428  for (j=k+1; j<n; ++j) {
429  u = kn+j; a[u] *= alinv;
430  }
431 
432  for (i=0; i<k; ++i) {
433  in = i*n;
434  for (j=0; j<n; ++j)
435  if (j!=k) {
436  u = in+j; a[u] -= a[in+k]*a[kn+j];
437  } // end if (j!=k)
438  }
439 
440  for (i=k+1; i<n; ++i) {
441  in = i*n;
442  for (j=0; j<n; ++j)
443  if (j!=k) {
444  u = in+j; a[u] -= a[in+k]*a[kn+j];
445  } // end if (j!=k)
446  }
447 
448  for (i=0; i<k; ++i) {
449  u=i*n+k; a[u] *= -alinv;
450  }
451 
452  for (i=k+1; i<n; ++i) {
453  u=i*n+k; a[u] *= -alinv;
454  }
455 
456  } // end for (k=0; k<n; ++k)
457 }
458 
473  const INT n)
474 {
475  INT i, j, k, l, ll, u;
476  INT icol, irow;
477  REAL vmax, dum, pivinv, temp;
478 
479  INT *work = (INT *)fasp_mem_calloc(3*n,sizeof(INT));
480  INT *indxc = work, *indxr = work+n, *ipiv = work+2*n;
481 
482  // ipiv, indxr, and indxc are used for book-keeping on the pivoting.
483  for ( j=0; j<n; j++ ) ipiv[j] = 0;
484 
485  // This is the main loop over the columns to be reduced.
486  for ( i=0; i<n; i++ ) {
487 
488  // This is the outer loop of the search for a pivot element.
489  vmax = 0.0;
490  for ( j=0; j<n; j++ ) {
491  if ( ipiv[j] != 1 ) {
492  for ( k=0; k<n; k++ ) {
493  if ( ipiv[k] == 0 ) {
494  u = j*n+k;
495  if ( ABS(a[u]) >= vmax ) {
496  vmax = ABS(a[u]); irow = j; icol = k;
497  }
498  }
499  } // end for k
500  }
501  } // end for j
502 
503  ++(ipiv[icol]);
504 
505  // We now have the pivot element, so we interchange rows, if needed, to put
506  // the pivot element on the diagonal. The columns are not physically
507  // interchanged, only relabeled: indxc[i], the column of the ith pivot
508  // element, is the ith column that is reduced, while indxr[i] is the row in
509  // which that pivot element was originally located. If indxr[i] ΜΈ= indxc[i]
510  // there is an implied column interchange. With this form of bookkeeping,
511  // the inverse matrix will be scrambled by columns.
512  if ( irow != icol ) {
513  for ( l=0; l<n; l++ ) SWAP(a[irow*n+l],a[icol*n+l]);
514  }
515 
516  indxr[i] = irow; indxc[i] = icol;
517  u = icol*n+icol;
518  if ( ABS(a[u]) < SMALLREAL ) {
519  printf("### WARNING: The matrix is nearly singular!\n");
520  }
521  pivinv = 1.0/a[u]; a[u]=1.0;
522  for ( l=0; l<n; l++ ) a[icol*n+l] *= pivinv;
523 
524  for ( ll=0; ll<n; ll++ ) {
525  if ( ll != icol ) {
526  u = ll*n+icol;
527  dum = a[u]; a[u] = 0.0;
528  for ( l=0; l<n; l++ ) a[ll*n+l] -= a[icol*n+l]*dum;
529  }
530  }
531  }
532  // This is the end of the main loop over columns of the reduction.
533 
534  // It only remains to unscramble the matrix in view of the column interchanges.
535  for ( l=n-1; l>=0; l-- ) {
536  if ( indxr[l] != indxc[l] )
537  for ( k=0; k<n; k++ ) SWAP(a[k*n+indxr[l]],a[k*n+indxc[l]]);
538  } // And we are done.
539 
540  fasp_mem_free(work);
541 }
542 
555  const INT n)
556 {
557  switch (n) {
558 
559  case 2:
561  break;
562 
563  case 3:
565  break;
566 
567  case 4:
569  break;
570 
571  case 5:
573  break;
574 
575  default:
577  break;
578 
579  }
580 
581  return FASP_SUCCESS;
582 }
583 
596  const INT n)
597 {
598 
599  REAL norm = 0.0;
600  REAL value = 0.0;
601 
602  INT i,j;
603 
604  for (i=0; i<n; i++){
605 
606  value = 0.0;
607 
608  for (j=0; j<n; j++){
609  value = value + ABS(A[i*n+j]);
610  }
611 
612  norm = MAX(norm, value);
613  }
614 
615  return norm;
616 }
617 
629 {
630  INT i;
631 
632  if (A==NULL) return;
633 
634  for (i=0;i<A->row;++i) fasp_mem_free(A->val[i]);
635  fasp_mem_free(A->val); A->val = NULL; A->row = 0;
636 }
637 
649 {
650  memset(a, 0X0, 4*sizeof(REAL));
651 
652  a[0] = 1.0; a[3] = 1.0;
653 }
654 
666 {
667  memset(a, 0X0, 9*sizeof(REAL));
668 
669  a[0] = 1.0; a[4] = 1.0; a[8] = 1.0;
670 }
671 
683 {
684  memset(a, 0X0, 25*sizeof(REAL));
685 
686  a[0] = 1.0;
687  a[6] = 1.0;
688  a[12] = 1.0;
689  a[18] = 1.0;
690  a[24] = 1.0;
691 }
692 
704 {
705  memset(a, 0X0, 49*sizeof(REAL));
706 
707  a[0] = 1.0;
708  a[8] = 1.0;
709  a[16] = 1.0;
710  a[24] = 1.0;
711  a[32] = 1.0;
712  a[40] = 1.0;
713  a[48] = 1.0;
714 }
715 
729  const INT n,
730  const INT n2)
731 {
732  memset(a, 0X0, n2*sizeof(REAL));
733 
734  switch (n) {
735 
736  case 2: {
737  a[0] = 1.0;
738  a[3] = 1.0;
739  }
740  break;
741 
742  case 3: {
743  a[0] = 1.0;
744  a[4] = 1.0;
745  a[8] = 1.0;
746  }
747  break;
748 
749  case 4: {
750  a[0] = 1.0;
751  a[5] = 1.0;
752  a[10] = 1.0;
753  a[15] = 1.0;
754  }
755  break;
756 
757  case 5: {
758  a[0] = 1.0;
759  a[6] = 1.0;
760  a[12] = 1.0;
761  a[18] = 1.0;
762  a[24] = 1.0;
763  }
764  break;
765 
766  case 6: {
767  a[0] = 1.0;
768  a[7] = 1.0;
769  a[14] = 1.0;
770  a[21] = 1.0;
771  a[28] = 1.0;
772  a[35] = 1.0;
773  }
774  break;
775 
776  case 7: {
777  a[0] = 1.0;
778  a[8] = 1.0;
779  a[16] = 1.0;
780  a[24] = 1.0;
781  a[32] = 1.0;
782  a[40] = 1.0;
783  a[48] = 1.0;
784  }
785  break;
786 
787  default: {
788  INT l;
789  for (l = 0; l < n; l ++) a[l*n+l] = 1.0;
790  }
791  break;
792  }
793 
794 }
795 
796 /*---------------------------------*/
797 /*-- End of File --*/
798 /*---------------------------------*/
void fasp_smat_identity_nc3(REAL *a)
Set a 3*3 full matrix to be a identity.
Definition: smat.c:665
#define REAL
Definition: fasp.h:67
void fasp_smat_identity_nc5(REAL *a)
Set a 5*5 full matrix to be a identity.
Definition: smat.c:682
INT fasp_blas_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: smat.c:554
void fasp_blas_smat_invp_nc(REAL *a, const INT n)
Compute the inverse of a matrix using Gauss Elimination with Pivoting.
Definition: smat.c:472
INT ** val
actual matrix entries
Definition: fasp.h:136
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
#define INT
Definition: fasp.h:64
INT row
number of rows
Definition: fasp.h:130
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_iden_free(idenmat *A)
Free idenmat sparse matrix data memeory space.
Definition: smat.c:628
Main header file for FASP.
void fasp_blas_smat_inv_nc7(REAL *a)
Compute the inverse matrix of a 7*7 matrix a.
Definition: smat.c:389
void fasp_blas_smat_inv_nc5(REAL *a)
Compute the inverse matrix of a 5*5 full matrix A (in place)
Definition: smat.c:173
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:55
#define SWAP(a, b)
Definition: smat.c:9
void fasp_blas_smat_inv_nc2(REAL *a)
Compute the inverse matrix of a 2*2 full matrix A (in place)
Definition: smat.c:25
REAL fasp_blas_smat_Linfinity(REAL *A, const INT n)
Compute the L infinity norm of A.
Definition: smat.c:595
void fasp_smat_identity(REAL *a, const INT n, const INT n2)
Set a n*n full matrix to be a identity.
Definition: smat.c:728
void fasp_blas_smat_inv_nc4(REAL *a)
Compute the inverse matrix of a 4*4 full matrix A (in place)
Definition: smat.c:115
void fasp_blas_smat_inv_nc3(REAL *a)
Compute the inverse matrix of a 3*3 full matrix A (in place)
Definition: smat.c:61
#define ABS(a)
Definition: fasp.h:74
void fasp_smat_identity_nc2(REAL *a)
Set a 2*2 full matrix to be a identity.
Definition: smat.c:648
void fasp_blas_smat_inv_nc(REAL *a, const INT n)
Compute the inverse of a matrix using Gauss Elimination.
Definition: smat.c:405
Dense matrix of INT type.
Definition: fasp.h:127
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define SMALLREAL
Definition: fasp_const.h:238
void fasp_smat_identity_nc7(REAL *a)
Set a 7*7 full matrix to be a identity.
Definition: smat.c:703