/[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 682 - (hide annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years, 7 months ago) by robwdcock
Original Path: trunk/paso/src/Solvers/Solver_SchurComplement.c
File MIME type: text/plain
File size: 10864 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26