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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 7 months 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 jgs 150 /* $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