"Fossies" - the Fresh Open Source Software Archive 
Member "laspack/precond.c" (27 Mar 1995, 3603 Bytes) of package /linux/privat/old/laspack.tgz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
1 /****************************************************************************/
2 /* precond.c */
3 /****************************************************************************/
4 /* */
5 /* PRECONDitioners for iterative solvers of systems of linear equations */
6 /* */
7 /* Copyright (C) 1992-1995 Tomas Skalicky. All rights reserved. */
8 /* */
9 /****************************************************************************/
10 /* */
11 /* ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS */
12 /* OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H) */
13 /* */
14 /****************************************************************************/
15
16 #include <math.h>
17
18 #include "laspack/precond.h"
19 #include "laspack/errhandl.h"
20 #include "laspack/qmatrix.h"
21 #include "laspack/operats.h"
22 #include "laspack/factor.h"
23 #include "laspack/itersolv.h"
24 #include "laspack/copyrght.h"
25
26 Vector *JacobiPrecond(QMatrix *A, Vector *y, Vector *c, double Omega)
27 /* Jacobi preconditioner */
28 {
29 Q_Lock(A);
30 V_Lock(y);
31 V_Lock(c);
32
33 /*
34 * Diag(A) / Omega y = c
35 *
36 * y = Omega Diag(A)^(-1) c
37 */
38
39 Asgn_VV(y, Mul_SV(Omega, MulInv_QV(Diag_Q(A), c)));
40
41 Q_Unlock(A);
42 V_Unlock(y);
43 V_Unlock(c);
44
45 return(y);
46 }
47
48 Vector *SSORPrecond(QMatrix *A, Vector *y, Vector *c, double Omega)
49 /* SSOR preconditioner */
50 {
51 Q_Lock(A);
52 V_Lock(y);
53 V_Lock(c);
54
55 /*
56 * 1 / (2 - Omega) * (Diag(A) / Omega + Lower(A)) * (Diag(A) / Omega)^(-1)
57 * * (Diag(A) / Omega + Upper(A)) y = c
58 *
59 * y = (2 - Omega) / Omega * (Diag(A) / Omega + Upper(A))^(-1) * Diag(A)
60 * * (Diag(A) / Omega + Lower(A))^(-1) c
61 */
62
63 Asgn_VV(y, Mul_SV((2.0 - Omega) / Omega,
64 MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Upper_Q(A)),
65 Mul_QV(Diag_Q(A),
66 MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A)), c)))));
67
68 Q_Unlock(A);
69 V_Unlock(y);
70 V_Unlock(c);
71
72 return(y);
73 }
74
75 Vector *ILUPrecond(QMatrix *A, Vector *y, Vector *c, double Omega)
76 /* incomplete factorization preconditioner */
77 {
78 Q_Lock(A);
79 V_Lock(y);
80 V_Lock(c);
81
82 /*
83 * if Q symmetric with rowwise stored elements
84 *
85 * (Diag(C) + Upper(C)) * Diag(C)^(-1) * (Diag(C) + Lower(C)) y = c
86 * y = (Diag(C) + Upper(C))^(-1) * Diag(C) * (Diag(C) + Lower(C))^(-1) c
87 *
88 * otherwise
89 *
90 * (Diag(B) + Lower(B)) * Diag(B)^(-1) * (Diag(B) + Upper(B)) y = c
91 * y = (Diag(B) + Upper(B))^(-1) * Diag(B) * (Diag(B) + Lower(B))^(-1) c
92 *
93 * B denotes the incomplete factorized matrix A
94 */
95
96 if (Q_GetSymmetry(A) && Q_GetElOrder(A) == Rowws)
97 Asgn_VV(y, MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Lower_Q(ILUFactor(A))),
98 Mul_QV(Diag_Q(ILUFactor(A)),
99 MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Upper_Q(ILUFactor(A))), c))));
100 else
101 Asgn_VV(y, MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Upper_Q(ILUFactor(A))),
102 Mul_QV(Diag_Q(ILUFactor(A)),
103 MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Lower_Q(ILUFactor(A))), c))));
104
105 Q_Unlock(A);
106 V_Unlock(y);
107 V_Unlock(c);
108
109 return(y);
110 }