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

Contents of /trunk/paso/src/SystemMatrix_MatrixVector.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years ago) by jgs
File MIME type: text/plain
File size: 6638 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Paso: matrix vector product with sparse matrix */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003,2004,2005 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Paso.h"
15 #include "SystemMatrix.h"
16
17 /**************************************************************************/
18
19 /* raw scaled vector update operation: out = alpha * A * in + beta * out */
20
21 /* has to be called within a parallel region */
22 /* barrier synconization is performed to make sure that the input vector available */
23
24 void Paso_SystemMatrix_MatrixVector(double alpha,
25 Paso_SystemMatrix* A,
26 double* in,
27 double beta,
28 double* out) {
29
30 register index_t ir,icol,iptr,icb,irb,irow,ic;
31 register double reg,reg1,reg2,reg3;
32 #pragma omp barrier
33
34 if (ABS(beta)>0.) {
35 #pragma omp for private(irow) schedule(static)
36 for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
37 out[irow] *= beta;
38 } else {
39 #pragma omp for private(irow) schedule(static)
40 for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
41 out[irow] = 0;
42 }
43
44 /* do the operation: */
45 if (ABS(alpha)>0) {
46 if (A ->col_block_size==1 && A->row_block_size ==1) {
47 switch(A->type) {
48 case CSR:
49 #pragma omp for private(irow,iptr,reg) schedule(static)
50 for (irow=0;irow< A->pattern->n_ptr;++irow) {
51 reg=0.;
52 for (iptr=(A->pattern->ptr[irow])-PTR_OFFSET;iptr<(A->pattern->ptr[irow+1])-PTR_OFFSET; ++iptr) {
53 reg += A->val[iptr] * in[A->pattern->index[iptr]-INDEX_OFFSET];
54 }
55 out[irow] += alpha * reg;
56 }
57 break;
58 case CSC:
59 /* TODO: parallelize (good luck!) */
60 #pragma omp single
61 for (icol=0;icol< A->pattern->n_ptr;++icol) {
62 for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; ++iptr) {
63 out[A->pattern->index[iptr]-INDEX_OFFSET]+= alpha * A->val[iptr] * in[icol];
64 }
65 }
66 break;
67 default:
68 Paso_setError(TYPE_ERROR,"Unknown matrix type in MVM.");
69 } /* switch A->type */
70 } else if (A ->col_block_size==2 && A->row_block_size ==2) {
71 switch(A->type) {
72 case CSR:
73 #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2) schedule(static)
74 for (ir=0;ir< A->pattern->n_ptr;ir++) {
75 reg1=0.;
76 reg2=0.;
77 for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
78 ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
79 reg1 += A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic];
80 reg2 += A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic];
81 }
82 out[ 2*ir] += alpha * reg1;
83 out[1+2*ir] += alpha * reg2;
84 }
85 break;
86 case CSC:
87 /* TODO: parallelize */
88 #pragma omp single
89 for (ic=0;ic< A->pattern->n_ptr;ic++) {
90 for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
91 ic=2*(A->pattern->index[iptr]-INDEX_OFFSET);
92 out[ 2*ir] += alpha * ( A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic] );
93 out[1+2*ir] += alpha * ( A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic] );
94 }
95 }
96 default:
97 Paso_setError(TYPE_ERROR,"Unknown matrix type in MVM.");
98 } /* switch A->type */
99 } else if (A ->col_block_size==3 && A->row_block_size ==3) {
100 switch(A->type) {
101 case CSR:
102 #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg1,reg2,reg3) schedule(static)
103 for (ir=0;ir< A->pattern->n_ptr;ir++) {
104 reg1=0.;
105 reg2=0.;
106 reg3=0.;
107 for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
108 ic=3*(A->pattern->index[iptr]-INDEX_OFFSET);
109 reg1 += A->val[iptr*9 ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic];
110 reg2 += A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic];
111 reg3 += A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic];
112 }
113 out[ 3*ir] += alpha * reg1;
114 out[1+3*ir] += alpha * reg2;
115 out[2+3*ir] += alpha * reg3;
116 }
117 break;
118 case CSC:
119 /* TODO: parallelize */
120 #pragma omp single
121 for (ic=0;ic< A->pattern->n_ptr;ic++) {
122 for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
123 ir=3*(A->pattern->index[iptr]-INDEX_OFFSET);
124 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] );
125 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] );
126 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] );
127 }
128 }
129 default:
130 Paso_setError(TYPE_ERROR,"Unknown matrix type in MVM.");
131 } /* switch A->type */
132 } else {
133 switch(A->type) {
134 case CSR:
135 #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)
136 for (ir=0;ir< A->pattern->n_ptr;ir++) {
137 for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
138 for (irb=0;irb< A->row_block_size;irb++) {
139 irow=irb+A->row_block_size*ir;
140 reg=0.;
141 for (icb=0;icb< A->col_block_size;icb++) {
142 icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
143 reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
144 }
145 out[irow] += alpha * reg;
146 }
147 }
148 }
149 break;
150 case CSC:
151 /* TODO: parallelize */
152 #pragma omp single
153 for (ic=0;ic< A->pattern->n_ptr;ic++) {
154 for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
155 for (irb=0;irb< A->row_block_size;irb++) {
156 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
157 for (icb=0;icb< A->col_block_size;icb++) {
158 icol=icb+A->col_block_size*ic;
159 out[irow] += alpha * A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
160 }
161 }
162 }
163 }
164 default:
165 Paso_setError(TYPE_ERROR,"Unknown matrix type in MVM.");
166 } /* switch A->type */
167 }
168 }
169 return;
170 }
171
172 /*
173 * $Log$
174 * Revision 1.2 2005/09/15 03:44:39 jgs
175 * Merge of development branch dev-02 back to main trunk on 2005-09-15
176 *
177 * Revision 1.1.2.1 2005/09/05 06:29:47 gross
178 * These files have been extracted from finley to define a stand alone libray for iterative
179 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
180 * has not been tested yet.
181 *
182 *
183 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26