/[escript]/trunk/paso/src/Solver_SchurComplement.c
ViewVC logotype

Annotation of /trunk/paso/src/Solver_SchurComplement.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/paso/src/Solvers/Solver_SchurComplement.c
File MIME type: text/plain
File size: 10213 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 150 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Paso: updates A_CC <- ACC-ACF AFF^{-1} AFC */
6    
7     /* no check of consistency of matrices !!!! */
8    
9     /**************************************************************/
10    
11     /* Copyrights by ACcESS Australia 2003,2004,2005 */
12     /* Author: gross@access.edu.au */
13    
14     /**************************************************************/
15    
16     #include "Paso.h"
17     #include "SystemMatrix.h"
18     #include "Solver.h"
19    
20     /**************************************************************/
21    
22    
23    
24     void Paso_Solver_updateIncompleteSchurComplement(Paso_SystemMatrix* A_CC,Paso_SystemMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot,Paso_SystemMatrix *A_FC) {
25     index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC;
26     dim_t index_CC_len;
27     bool_t set_A;
28     dim_t n_rows=A_CC->num_rows;
29     dim_t n_block=A_CC->row_block_size;
30     register double A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,
31     invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,
32     A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33;
33     if (n_block==1) {
34     #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,col_CF,set_A,iPtr_CF,iPtr_FC,col_FC,where_p,A11) schedule(static)
35     for (i = 0; i < n_rows;++i) {
36     iPtr_CC=A_CC->pattern->ptr[i];
37     index_CC=&(A_CC->pattern->index[iPtr_CC]);
38     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
39     /* now we run through the columns of A_CF in row i */
40     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
41     col_CF=A_CF->pattern->index[iPtr_CF];
42     set_A=TRUE;
43     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
44     col_FC=A_FC->pattern->index[iPtr_FC];
45     /* is (i,col_FC) in the shape of A_CC ? */
46     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
47     if (where_p!=NULL) {
48     if (set_A) {
49     A11=A_CF->val[iPtr_CF]*invA_FF[col_CF];
50     set_A=FALSE;
51     }
52     A_CC->val[iPtr_CC+(index_t)(where_p-index_CC)]-=A11*A_FC->val[iPtr_FC];
53     }
54     } /* end of iPtr_FC loop */
55     } /* end of iPtr_CF loop */
56     } /* end of irow loop */
57     } else if (n_block==2) {
58     #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_12,A_CF_22,invA_FF_11,invA_FF_21,invA_FF_12,invA_FF_22,A11,A21,A12,A22,A_FC_11,A_FC_21,A_FC_12,A_FC_22) schedule(static)
59     for (i = 0; i < n_rows;++i) {
60     iPtr_CC=A_CC->pattern->ptr[i];
61     index_CC=&(A_CC->pattern->index[iPtr_CC]);
62     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
63     /* now we run through the columns of A_CF in row i */
64     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
65     col_CF=A_CF->pattern->index[iPtr_CF];
66     set_A=TRUE;
67     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
68     col_FC=A_FC->pattern->index[iPtr_FC];
69     /* is (i,col_FC) in the shape of A_CC ? */
70     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
71     if (where_p!=NULL) {
72     iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC);
73     /* this calculutes A_CF*invA_FF(i,col_CF) */
74     if (set_A) {
75     A_CF_11=A_CF->val[4*iPtr_CF ];
76     A_CF_21=A_CF->val[4*iPtr_CF+1];
77     A_CF_12=A_CF->val[4*iPtr_CF+2];
78     A_CF_22=A_CF->val[4*iPtr_CF+3];
79    
80     invA_FF_11=invA_FF[4*col_CF ];
81     invA_FF_21=invA_FF[4*col_CF+1];
82     invA_FF_12=invA_FF[4*col_CF+2];
83     invA_FF_22=invA_FF[4*col_CF+3];
84    
85     A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21;
86     A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21;
87     A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22;
88     A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22;
89    
90     set_A=FALSE;
91     }
92    
93     A_FC_11=A_FC->val[4*iPtr_FC ];
94     A_FC_21=A_FC->val[4*iPtr_FC+1];
95     A_FC_12=A_FC->val[4*iPtr_FC+2];
96     A_FC_22=A_FC->val[4*iPtr_FC+3];
97    
98     A_CC->val[4*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21;
99     A_CC->val[4*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21;
100     A_CC->val[4*iPtr_CC_2+2]-=A11*A_FC_12+A12*A_FC_22;
101     A_CC->val[4*iPtr_CC_2+3]-=A21*A_FC_12+A22*A_FC_22;
102    
103     }
104     } /* end of iPtr_FC loop */
105     } /* end of iPtr_CF loop */
106     } /* end of irow loop */
107     } else if (n_block==3) {
108     #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33) schedule(static)
109     for (i = 0; i < n_rows;++i) {
110     iPtr_CC=A_CC->pattern->ptr[i];
111     index_CC=&(A_CC->pattern->index[iPtr_CC]);
112     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
113     /* now we run through the columns of A_CF in row i */
114     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
115     col_CF=A_CF->pattern->index[iPtr_CF];
116     set_A=TRUE;
117     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
118     col_FC=A_FC->pattern->index[iPtr_FC];
119     /* is (i,col_FC) in the shape of A_CC ? */
120     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
121     if (where_p!=NULL) {
122     iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC);
123     /* this calculutes A_CF*invA_FF(i,col_CF) */
124     if (set_A) {
125     A_CF_11=A_CF->val[9*iPtr_CF ];
126     A_CF_21=A_CF->val[9*iPtr_CF+1];
127     A_CF_31=A_CF->val[9*iPtr_CF+2];
128     A_CF_12=A_CF->val[9*iPtr_CF+3];
129     A_CF_22=A_CF->val[9*iPtr_CF+4];
130     A_CF_32=A_CF->val[9*iPtr_CF+5];
131     A_CF_13=A_CF->val[9*iPtr_CF+6];
132     A_CF_23=A_CF->val[9*iPtr_CF+7];
133     A_CF_33=A_CF->val[9*iPtr_CF+8];
134    
135     invA_FF_11=invA_FF[9*col_CF ];
136     invA_FF_21=invA_FF[9*col_CF+1];
137     invA_FF_31=invA_FF[9*col_CF+2];
138     invA_FF_12=invA_FF[9*col_CF+3];
139     invA_FF_22=invA_FF[9*col_CF+4];
140     invA_FF_32=invA_FF[9*col_CF+5];
141     invA_FF_13=invA_FF[9*col_CF+6];
142     invA_FF_23=invA_FF[9*col_CF+7];
143     invA_FF_33=invA_FF[9*col_CF+8];
144    
145     A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21+A_CF_13*invA_FF_31;
146     A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21+A_CF_23*invA_FF_31;
147     A31=A_CF_31*invA_FF_11+A_CF_32*invA_FF_21+A_CF_33*invA_FF_31;
148     A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22+A_CF_13*invA_FF_32;
149     A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22+A_CF_23*invA_FF_32;
150     A32=A_CF_31*invA_FF_12+A_CF_32*invA_FF_22+A_CF_33*invA_FF_32;
151     A13=A_CF_11*invA_FF_13+A_CF_12*invA_FF_23+A_CF_13*invA_FF_33;
152     A23=A_CF_21*invA_FF_13+A_CF_22*invA_FF_23+A_CF_23*invA_FF_33;
153     A33=A_CF_31*invA_FF_13+A_CF_32*invA_FF_23+A_CF_33*invA_FF_33;
154    
155     set_A=FALSE;
156     }
157    
158     A_FC_11=A_FC->val[9*iPtr_FC ];
159     A_FC_21=A_FC->val[9*iPtr_FC+1];
160     A_FC_31=A_FC->val[9*iPtr_FC+2];
161     A_FC_12=A_FC->val[9*iPtr_FC+3];
162     A_FC_22=A_FC->val[9*iPtr_FC+4];
163     A_FC_32=A_FC->val[9*iPtr_FC+5];
164     A_FC_13=A_FC->val[9*iPtr_FC+6];
165     A_FC_23=A_FC->val[9*iPtr_FC+7];
166     A_FC_33=A_FC->val[9*iPtr_FC+8];
167    
168     A_CC->val[9*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21+A13*A_FC_31;
169     A_CC->val[9*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21+A23*A_FC_31;
170     A_CC->val[9*iPtr_CC_2+2]-=A31*A_FC_11+A32*A_FC_21+A33*A_FC_31;
171     A_CC->val[9*iPtr_CC_2+3]-=A11*A_FC_12+A12*A_FC_22+A13*A_FC_32;
172     A_CC->val[9*iPtr_CC_2+4]-=A21*A_FC_12+A22*A_FC_22+A23*A_FC_32;
173     A_CC->val[9*iPtr_CC_2+5]-=A31*A_FC_12+A32*A_FC_22+A33*A_FC_32;
174     A_CC->val[9*iPtr_CC_2+6]-=A11*A_FC_13+A12*A_FC_23+A13*A_FC_33;
175     A_CC->val[9*iPtr_CC_2+7]-=A21*A_FC_13+A22*A_FC_23+A23*A_FC_33;
176     A_CC->val[9*iPtr_CC_2+8]-=A31*A_FC_13+A32*A_FC_23+A33*A_FC_33;
177    
178     }
179     } /* end of iPtr_FC loop */
180     } /* end of iPtr_CF loop */
181     } /* end of irow loop */
182     }
183     }
184     /*
185     * $Log$
186     * Revision 1.2 2005/09/15 03:44:40 jgs
187     * Merge of development branch dev-02 back to main trunk on 2005-09-15
188     *
189     * Revision 1.1.2.1 2005/09/05 06:29:50 gross
190     * These files have been extracted from finley to define a stand alone libray for iterative
191     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
192     * has not been tested yet.
193     *
194     *
195     */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26