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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 115 - (show annotations)
Fri Mar 4 07:12:47 2005 UTC (14 years, 8 months ago) by jgs
File MIME type: text/plain
File size: 7944 byte(s)
*** empty log message ***

1 /* $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 register maybelong ir,icol,iptr,icb,irb,irow,ic;
74 register double reg,reg1,reg2,reg3;
75 #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 #pragma omp for private(irow,iptr,reg) schedule(static)
93 for (irow=0;irow< A->pattern->n_ptr;++irow) {
94 reg=0.;
95 for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {
96 reg += A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];
97 }
98 out[irow] += alpha * reg;
99 }
100 break;
101 case CSC:
102 /* TODO: parallelize (good luck!) */
103 #pragma omp single
104 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 }
108 }
109 break;
110 default:
111 Finley_ErrorCode=TYPE_ERROR;
112 sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
113 } /* switch A->type */
114 } 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 } else {
179 switch(A->type) {
180 case CSR:
181 #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
182 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 for (irb=0;irb< A->row_block_size;irb++) {
185 irow=irb+A->row_block_size*ir;
186 reg=0.;
187 for (icb=0;icb< A->col_block_size;icb++) {
188 icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
189 reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
190 }
191 out[irow] += alpha * reg;
192 }
193 }
194 }
195 break;
196 case CSC:
197 /* TODO: parallelize */
198 #pragma omp single
199 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 for (irb=0;irb< A->row_block_size;irb++) {
202 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
203 for (icb=0;icb< A->col_block_size;icb++) {
204 icol=icb+A->col_block_size*ic;
205 out[irow] += alpha * A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
206 }
207 }
208 }
209 }
210 default:
211 Finley_ErrorCode=TYPE_ERROR;
212 sprintf(Finley_ErrorMsg,"Unknown matrix type in MVM.");
213 } /* switch A->type */
214 }
215 }
216 return;
217 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26