/[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 100 - (show annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 5610 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 } else if (A ->symmetric) {
50 Finley_ErrorCode=SYSTEM_ERROR;
51 sprintf(Finley_ErrorMsg,"matrix-vector product for symmetrix matric has not been implemented yet.");
52 }
53
54 /* delegate to routine */
55
56 if (Finley_ErrorCode==NO_ERROR) {
57 #pragma omp parallel
58 Finley_RawScaledSystemMatrixVector(alpha, A, getSampleData(in,0), beta, getSampleData(out,0));
59 }
60 return;
61 }
62
63 /**************************************************************************/
64
65 /* raw scaled vector update operation: out = alpha * A * in + beta * out */
66
67 /* has to be called within a parallel region */
68 /* barrier synconization is performed to make sure that the input vector available */
69
70 void Finley_RawScaledSystemMatrixVector(double alpha,
71 Finley_SystemMatrix* A,
72 double* in,
73 double beta,
74 double* out) {
75
76 maybelong ir,icol,iptr,icb,irb,irow,ic;
77 #pragma omp barrier
78
79 if (ABS(beta)>0.) {
80 #pragma omp for private(irow) schedule(static)
81 for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
82 out[irow] *= beta;
83 } else {
84 #pragma omp for private(irow) schedule(static)
85 for (irow=0;irow < A->num_rows * A->row_block_size;irow++)
86 out[irow] = 0;
87 }
88
89 /* do the operation: */
90 if (ABS(alpha)>0) {
91 if (A ->col_block_size==1 && A->row_block_size ==1) {
92 switch(A->type) {
93 case CSR:
94 #pragma omp for private(irow,iptr) schedule(static)
95 for (irow=0;irow< A->num_rows;irow++) {
96 for (iptr=A->ptr[irow]-PTR_OFFSET;iptr<A->ptr[irow+1]-PTR_OFFSET; iptr++) {
97 out[irow] += alpha * A->val[iptr] * in[A->index[iptr]-INDEX_OFFSET];
98 }
99 }
100 break;
101 case CSC:
102 /* TODO: parallelize */
103 #pragma omp single
104 for (icol=0;icol< A->num_cols;icol++) {
105 for (iptr=A->ptr[icol]-PTR_OFFSET;iptr<A->ptr[icol+1]-PTR_OFFSET; iptr++) {
106 out[A->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.");
113 } /* switch A->type */
114 } else {
115 switch(A->type) {
116 case CSR:
117 #pragma omp for private(ir,iptr,irb,icb,irow,icol) schedule(static)
118 for (ir=0;ir< A->num_rows;ir++) {
119 for (iptr=A->ptr[ir]-PTR_OFFSET;iptr<A->ptr[ir+1]-PTR_OFFSET; iptr++) {
120 for (irb=0;irb< A->row_block_size;irb++) {
121 irow=irb+A->row_block_size*ir;
122 for (icb=0;icb< A->col_block_size;icb++) {
123 icol=icb+A->col_block_size*(A->index[iptr]-INDEX_OFFSET);
124 out[irow] += alpha * A->val[iptr*A->row_block_size*A->col_block_size+icb+A->col_block_size*irb] * in[icol];
125 }
126 }
127 }
128 }
129 break;
130 case CSC:
131 /* TODO: parallelize */
132 #pragma omp single
133 for (ic=0;ic< A->num_cols;ic++) {
134 for (iptr=A->ptr[ic]-PTR_OFFSET;iptr<A->ptr[ic+1]-PTR_OFFSET; iptr++) {
135 for (irb=0;irb< A->row_block_size;irb++) {
136 irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);
137 for (icb=0;icb< A->col_block_size;icb++) {
138 icol=icb+A->col_block_size*ic;
139 out[irow] += alpha * A->val[iptr*A->row_block_size*A->col_block_size+icb+A->col_block_size*irb] * in[icol];
140 }
141 }
142 }
143 }
144 default:
145 Finley_ErrorCode=TYPE_ERROR;
146 sprintf(Finley_ErrorMsg,"Unknown matrix type.");
147 } /* switch A->type */
148 }
149 }
150 return;
151 }
152
153 /*
154 * $Log$
155 * Revision 1.3 2004/12/15 03:48:46 jgs
156 * *** empty log message ***
157 *
158 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
159 * initial import of project esys2
160 *
161 * Revision 1.2.2.2 2004/10/26 06:36:39 jgs
162 * committing Lutz's changes to branch jgs
163 *
164 * Revision 1.3 2004/10/21 04:55:54 gross
165 * bug in CSC MVM fixed
166 *
167 * Revision 1.2 2004/08/13 00:12:53 gross
168 * Gradtest is getting further now. PDE assemblage has been added but not tested.
169 *
170 * Revision 1.1 2004/07/02 04:21:13 gross
171 * Finley C code has been included
172 *
173 *
174 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26