/[escript]/trunk/esys2/finley/src/finleyC/System_MVM.c
ViewVC logotype

Annotation of /trunk/esys2/finley/src/finleyC/System_MVM.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 9198 byte(s)
Merge of development branch back to main trunk on 2005-07-08

1 jgs 82 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Finley: matrix vector product with sparse matrix */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2003/04 */
10     /* Author: gross@access.edu.au */
11    
12     /**************************************************************/
13    
14     #include "Finley.h"
15     #include "System.h"
16    
17     /**************************************************************/
18    
19     /* vector update operation: out=out+A*in */
20    
21     void Finley_SystemMatrixVector(escriptDataC * out,Finley_SystemMatrix* A, escriptDataC * in) {
22     Finley_ScaledSystemMatrixVector(1.0, A, in, 1.0, out);
23     }
24    
25     /**************************************************************/
26    
27     /* scaled vector update operation: out = alpha * A * in + beta * out */
28    
29     void Finley_ScaledSystemMatrixVector(double alpha,
30     Finley_SystemMatrix* A,
31     escriptDataC * in,
32     double beta,
33     escriptDataC * out) {
34     Finley_ErrorCode=NO_ERROR;
35     /* check structure: */
36    
37     if (!isExpanded(in)) {
38     Finley_ErrorCode=TYPE_ERROR;
39     sprintf(Finley_ErrorMsg,"matrix vector product : input Data object has to be expanded");
40     } else if (!isExpanded(out)) {
41     Finley_ErrorCode=TYPE_ERROR;
42     sprintf(Finley_ErrorMsg,"matrix vector product : output Data object has to be expanded");
43     } else if ( getLength(in) != A ->col_block_size * A ->num_cols) {
44     Finley_ErrorCode=TYPE_ERROR;
45     sprintf(Finley_ErrorMsg,"matrix vector product : lnegth of input Data object and matrix column dimensions don't match.");
46     } else if (getLength(out) != A ->row_block_size * A ->num_rows) {
47     Finley_ErrorCode=TYPE_ERROR;
48     sprintf(Finley_ErrorMsg,"matrix vector product : lnegth of output Data object and matrix row dimensions don't match.");
49     }
50    
51     /* delegate to routine */
52    
53     if (Finley_ErrorCode==NO_ERROR) {
54     #pragma omp parallel
55     Finley_RawScaledSystemMatrixVector(alpha, A, getSampleData(in,0), beta, getSampleData(out,0));
56     }
57     return;
58     }
59    
60     /**************************************************************************/
61    
62     /* raw scaled vector update operation: out = alpha * A * in + beta * out */
63    
64     /* has to be called within a parallel region */
65     /* barrier synconization is performed to make sure that the input vector available */
66    
67     void Finley_RawScaledSystemMatrixVector(double alpha,
68     Finley_SystemMatrix* A,
69     double* in,
70     double beta,
71     double* out) {
72    
73 jgs 123 register index_t ir,icol,iptr,icb,irb,irow,ic;
74 jgs 115 register double reg,reg1,reg2,reg3;
75 jgs 82 #pragma omp barrier
76    
77     if (ABS(beta)>0.) {
78     #pragma omp for private(irow) schedule(static)
79     for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
80     out[irow] *= beta;
81     } else {
82     #pragma omp for private(irow) schedule(static)
83     for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
84     out[irow] = 0;
85     }
86    
87     /* do the operation: */
88     if (ABS(alpha)>0) {
89     if (A ->col_block_size==1 && A->row_block_size ==1) {
90     switch(A->type) {
91     case CSR:
92 jgs 113 #pragma omp for private(irow,iptr,reg) schedule(static)
93 jgs 102 for (irow=0;irow< A->pattern->n_ptr;++irow) {
94 jgs 113 reg=0.;
95 jgs 102 for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {
96 jgs 113 reg += A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];
97 jgs 82 }
98 jgs 113 out[irow] += alpha * reg;
99 jgs 82 }
100     break;
101     case CSC:
102 jgs 115 /* TODO: parallelize (good luck!) */
103 jgs 82 #pragma omp single
104 jgs 102 for (icol=0;icol< A->pattern->n_ptr;++icol) {
105     for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; ++iptr) {
106     out[A->pattern->index[iptr]-INDEX_OFFSET]+= alpha * A->val[iptr] * in[icol];
107 jgs 82 }
108     }
109     break;
110     default:
111     Finley_ErrorCode=TYPE_ERROR;
112 jgs 102 sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
113 jgs 82 } /* switch A->type */
114 jgs 113 } else if (A ->col_block_size==2 && A->row_block_size ==2) {
115     switch(A->type) {
116     case CSR:
117     #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2) schedule(static)
118     for (ir=0;ir< A->pattern->n_ptr;ir++) {
119     reg1=0.;
120     reg2=0.;
121     for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
122     ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
123     reg1 += A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic];
124     reg2 += A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic];
125     }
126     out[ 2*ir] += alpha * reg1;
127     out[1+2*ir] += alpha * reg2;
128     }
129     break;
130     case CSC:
131     /* TODO: parallelize */
132     #pragma omp single
133     for (ic=0;ic< A->pattern->n_ptr;ic++) {
134     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
135     ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
136     out[ 2*ir] += alpha * ( A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );
137     out[1+2*ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );
138     }
139     }
140     default:
141     Finley_ErrorCode=TYPE_ERROR;
142     sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
143     } /* switch A->type */
144     } else if (A ->col_block_size==3 && A->row_block_size ==3) {
145     switch(A->type) {
146     case CSR:
147     #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2,reg3) schedule(static)
148     for (ir=0;ir< A->pattern->n_ptr;ir++) {
149     reg1=0.;
150     reg2=0.;
151     reg3=0.;
152     for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
153     ic=3*(A->pattern->index[iptr]-INDEX_OFFSET);
154     reg1 += A->val[iptr*9 ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic];
155     reg2 += A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic];
156     reg3 += A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic];
157     }
158     out[ 3*ir] += alpha * reg1;
159     out[1+3*ir] += alpha * reg2;
160     out[2+3*ir] += alpha * reg3;
161     }
162     break;
163     case CSC:
164     /* TODO: parallelize */
165     #pragma omp single
166     for (ic=0;ic< A->pattern->n_ptr;ic++) {
167     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
168     ir=3*(A->pattern->index[iptr]-INDEX_OFFSET);
169     out[ 3*ir] += alpha * ( A->val[iptr*9 ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic] );
170     out[1+3*ir] += alpha * ( A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic] );
171     out[2+3*ir] += alpha * ( A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic] );
172     }
173     }
174     default:
175     Finley_ErrorCode=TYPE_ERROR;
176     sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
177     } /* switch A->type */
178 jgs 82 } else {
179     switch(A->type) {
180     case CSR:
181 jgs 113 #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
182 jgs 102 for (ir=0;ir< A->pattern->n_ptr;ir++) {
183     for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
184 jgs 82 for (irb=0;irb< A->row_block_size;irb++) {
185     irow=irb+A->row_block_size*ir;
186 jgs 113 reg=0.;
187 jgs 82 for (icb=0;icb< A->col_block_size;icb++) {
188 jgs 102 icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
189 jgs 113 reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
190 jgs 82 }
191 jgs 113 out[irow] += alpha * reg;
192 jgs 82 }
193     }
194     }
195     break;
196     case CSC:
197     /* TODO: parallelize */
198     #pragma omp single
199 jgs 102 for (ic=0;ic< A->pattern->n_ptr;ic++) {
200     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
201 jgs 82 for (irb=0;irb< A->row_block_size;irb++) {
202 jgs 102 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
203 jgs 82 for (icb=0;icb< A->col_block_size;icb++) {
204     icol=icb+A->col_block_size*ic;
205 jgs 102 out[irow] += alpha * A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
206 jgs 82 }
207     }
208     }
209     }
210     default:
211     Finley_ErrorCode=TYPE_ERROR;
212 jgs 102 sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
213 jgs 82 } /* switch A->type */
214     }
215     }
216     return;
217     }
218 jgs 123
219     /*
220     * $Log$
221     * Revision 1.7 2005/07/08 04:07:57 jgs
222     * Merge of development branch back to main trunk on 2005-07-08
223     *
224     * Revision 1.1.1.1.2.5 2005/06/29 02:34:56 gross
225     * some changes towards 64 integers in finley
226     *
227     * Revision 1.1.1.1.2.4 2005/03/02 23:35:06 gross
228     * reimplementation of the ILU in Finley. block size>1 still needs some testing
229     *
230     * Revision 1.1.1.1.2.3 2005/02/17 03:23:02 gross
231     * some performance improvements in MVM
232     *
233     * Revision 1.1.1.1.2.2 2004/12/07 10:12:05 gross
234     * GMRES added
235     *
236     * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
237     * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
238     *
239     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
240     * initial import of project esys2
241     *
242     * Revision 1.2.2.2 2004/10/26 06:36:39 jgs
243     * committing Lutz's changes to branch jgs
244     *
245     * Revision 1.3 2004/10/21 04:55:54 gross
246     * bug in CSC MVM fixed
247     *
248     * Revision 1.2 2004/08/13 00:12:53 gross
249     * Gradtest is getting further now. PDE assemblage has been added but not tested.
250     *
251     * Revision 1.1 2004/07/02 04:21:13 gross
252     * Finley C code has been included
253     *
254     *
255     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26