/[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 1387 - (hide annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/paso/src/Solver_SchurComplement.c
File MIME type: text/plain
File size: 10657 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
1 ksteube 1312
2 jgs 150 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 dhawcroft 631
16 jgs 150 /**************************************************************/
17    
18     /* Paso: updates A_CC <- ACC-ACF AFF^{-1} AFC */
19    
20     /* no check of consistency of matrices !!!! */
21    
22     /**************************************************************/
23    
24     /* Copyrights by ACcESS Australia 2003,2004,2005 */
25     /* Author: gross@access.edu.au */
26    
27     /**************************************************************/
28    
29 gross 700 #include "Paso.h"
30 ksteube 1312 #include "SparseMatrix.h"
31 jgs 150 #include "Solver.h"
32    
33     /**************************************************************/
34    
35    
36    
37 ksteube 1312 void Paso_Solver_updateIncompleteSchurComplement(Paso_SparseMatrix* A_CC,Paso_SparseMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot,Paso_SparseMatrix *A_FC) {
38    
39 jgs 150 index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC;
40     dim_t index_CC_len;
41     bool_t set_A;
42 ksteube 1312 dim_t n_loc_rows=A_CC->numRows;
43 jgs 150 dim_t n_block=A_CC->row_block_size;
44     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,
45     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,
46     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;
47     if (n_block==1) {
48     #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)
49 ksteube 1312 for (i = 0; i < n_loc_rows;++i) {
50 jgs 150 iPtr_CC=A_CC->pattern->ptr[i];
51     index_CC=&(A_CC->pattern->index[iPtr_CC]);
52     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
53     /* now we run through the columns of A_CF in row i */
54     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
55     col_CF=A_CF->pattern->index[iPtr_CF];
56     set_A=TRUE;
57     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
58     col_FC=A_FC->pattern->index[iPtr_FC];
59     /* is (i,col_FC) in the shape of A_CC ? */
60     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
61     if (where_p!=NULL) {
62     if (set_A) {
63     A11=A_CF->val[iPtr_CF]*invA_FF[col_CF];
64     set_A=FALSE;
65     }
66     A_CC->val[iPtr_CC+(index_t)(where_p-index_CC)]-=A11*A_FC->val[iPtr_FC];
67     }
68     } /* end of iPtr_FC loop */
69     } /* end of iPtr_CF loop */
70     } /* end of irow loop */
71     } else if (n_block==2) {
72     #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)
73 ksteube 1312 for (i = 0; i < n_loc_rows;++i) {
74 jgs 150 iPtr_CC=A_CC->pattern->ptr[i];
75     index_CC=&(A_CC->pattern->index[iPtr_CC]);
76     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
77     /* now we run through the columns of A_CF in row i */
78     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
79     col_CF=A_CF->pattern->index[iPtr_CF];
80     set_A=TRUE;
81     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
82     col_FC=A_FC->pattern->index[iPtr_FC];
83     /* is (i,col_FC) in the shape of A_CC ? */
84     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
85     if (where_p!=NULL) {
86     iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC);
87     /* this calculutes A_CF*invA_FF(i,col_CF) */
88     if (set_A) {
89     A_CF_11=A_CF->val[4*iPtr_CF ];
90     A_CF_21=A_CF->val[4*iPtr_CF+1];
91     A_CF_12=A_CF->val[4*iPtr_CF+2];
92     A_CF_22=A_CF->val[4*iPtr_CF+3];
93    
94     invA_FF_11=invA_FF[4*col_CF ];
95     invA_FF_21=invA_FF[4*col_CF+1];
96     invA_FF_12=invA_FF[4*col_CF+2];
97     invA_FF_22=invA_FF[4*col_CF+3];
98    
99     A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21;
100     A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21;
101     A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22;
102     A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22;
103    
104     set_A=FALSE;
105     }
106    
107     A_FC_11=A_FC->val[4*iPtr_FC ];
108     A_FC_21=A_FC->val[4*iPtr_FC+1];
109     A_FC_12=A_FC->val[4*iPtr_FC+2];
110     A_FC_22=A_FC->val[4*iPtr_FC+3];
111    
112     A_CC->val[4*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21;
113     A_CC->val[4*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21;
114     A_CC->val[4*iPtr_CC_2+2]-=A11*A_FC_12+A12*A_FC_22;
115     A_CC->val[4*iPtr_CC_2+3]-=A21*A_FC_12+A22*A_FC_22;
116    
117     }
118     } /* end of iPtr_FC loop */
119     } /* end of iPtr_CF loop */
120     } /* end of irow loop */
121     } else if (n_block==3) {
122     #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)
123 ksteube 1312 for (i = 0; i < n_loc_rows;++i) {
124 jgs 150 iPtr_CC=A_CC->pattern->ptr[i];
125     index_CC=&(A_CC->pattern->index[iPtr_CC]);
126     index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
127     /* now we run through the columns of A_CF in row i */
128     for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) {
129     col_CF=A_CF->pattern->index[iPtr_CF];
130     set_A=TRUE;
131     for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
132     col_FC=A_FC->pattern->index[iPtr_FC];
133     /* is (i,col_FC) in the shape of A_CC ? */
134     where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
135     if (where_p!=NULL) {
136     iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC);
137     /* this calculutes A_CF*invA_FF(i,col_CF) */
138     if (set_A) {
139     A_CF_11=A_CF->val[9*iPtr_CF ];
140     A_CF_21=A_CF->val[9*iPtr_CF+1];
141     A_CF_31=A_CF->val[9*iPtr_CF+2];
142     A_CF_12=A_CF->val[9*iPtr_CF+3];
143     A_CF_22=A_CF->val[9*iPtr_CF+4];
144     A_CF_32=A_CF->val[9*iPtr_CF+5];
145     A_CF_13=A_CF->val[9*iPtr_CF+6];
146     A_CF_23=A_CF->val[9*iPtr_CF+7];
147     A_CF_33=A_CF->val[9*iPtr_CF+8];
148    
149     invA_FF_11=invA_FF[9*col_CF ];
150     invA_FF_21=invA_FF[9*col_CF+1];
151     invA_FF_31=invA_FF[9*col_CF+2];
152     invA_FF_12=invA_FF[9*col_CF+3];
153     invA_FF_22=invA_FF[9*col_CF+4];
154     invA_FF_32=invA_FF[9*col_CF+5];
155     invA_FF_13=invA_FF[9*col_CF+6];
156     invA_FF_23=invA_FF[9*col_CF+7];
157     invA_FF_33=invA_FF[9*col_CF+8];
158    
159     A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21+A_CF_13*invA_FF_31;
160     A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21+A_CF_23*invA_FF_31;
161     A31=A_CF_31*invA_FF_11+A_CF_32*invA_FF_21+A_CF_33*invA_FF_31;
162     A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22+A_CF_13*invA_FF_32;
163     A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22+A_CF_23*invA_FF_32;
164     A32=A_CF_31*invA_FF_12+A_CF_32*invA_FF_22+A_CF_33*invA_FF_32;
165     A13=A_CF_11*invA_FF_13+A_CF_12*invA_FF_23+A_CF_13*invA_FF_33;
166     A23=A_CF_21*invA_FF_13+A_CF_22*invA_FF_23+A_CF_23*invA_FF_33;
167     A33=A_CF_31*invA_FF_13+A_CF_32*invA_FF_23+A_CF_33*invA_FF_33;
168    
169     set_A=FALSE;
170     }
171    
172     A_FC_11=A_FC->val[9*iPtr_FC ];
173     A_FC_21=A_FC->val[9*iPtr_FC+1];
174     A_FC_31=A_FC->val[9*iPtr_FC+2];
175     A_FC_12=A_FC->val[9*iPtr_FC+3];
176     A_FC_22=A_FC->val[9*iPtr_FC+4];
177     A_FC_32=A_FC->val[9*iPtr_FC+5];
178     A_FC_13=A_FC->val[9*iPtr_FC+6];
179     A_FC_23=A_FC->val[9*iPtr_FC+7];
180     A_FC_33=A_FC->val[9*iPtr_FC+8];
181    
182     A_CC->val[9*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21+A13*A_FC_31;
183     A_CC->val[9*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21+A23*A_FC_31;
184     A_CC->val[9*iPtr_CC_2+2]-=A31*A_FC_11+A32*A_FC_21+A33*A_FC_31;
185     A_CC->val[9*iPtr_CC_2+3]-=A11*A_FC_12+A12*A_FC_22+A13*A_FC_32;
186     A_CC->val[9*iPtr_CC_2+4]-=A21*A_FC_12+A22*A_FC_22+A23*A_FC_32;
187     A_CC->val[9*iPtr_CC_2+5]-=A31*A_FC_12+A32*A_FC_22+A33*A_FC_32;
188     A_CC->val[9*iPtr_CC_2+6]-=A11*A_FC_13+A12*A_FC_23+A13*A_FC_33;
189     A_CC->val[9*iPtr_CC_2+7]-=A21*A_FC_13+A22*A_FC_23+A23*A_FC_33;
190     A_CC->val[9*iPtr_CC_2+8]-=A31*A_FC_13+A32*A_FC_23+A33*A_FC_33;
191    
192     }
193     } /* end of iPtr_FC loop */
194     } /* end of iPtr_CF loop */
195     } /* end of irow loop */
196     }
197     }
198     /*
199     * $Log$
200     * Revision 1.2 2005/09/15 03:44:40 jgs
201     * Merge of development branch dev-02 back to main trunk on 2005-09-15
202     *
203     * Revision 1.1.2.1 2005/09/05 06:29:50 gross
204     * These files have been extracted from finley to define a stand alone libray for iterative
205     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
206     * has not been tested yet.
207     *
208     *
209     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26