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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (15 years, 5 months ago) by jgs
File MIME type: text/plain
File size: 5230 byte(s)
*** empty log message ***

1 jgs 82 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Finley: SystemMatrix */
6    
7     /* nullify rows and columns in the matrix */
8    
9     /* the rows and columns are marked by positive values in */
10     /* mask_row and mask_col. Values on the main diagonal */
11     /* which are marked to set to zero by both mask_row and */
12     /* mask_col are set to main_diagonal_value */
13    
14    
15     /**************************************************************/
16    
17     /* Copyrights by ACcESS Australia 2003 */
18     /* Author: gross@access.edu.au */
19    
20     /**************************************************************/
21    
22     #include "escript/Data/DataC.h"
23     #include "Finley.h"
24     #include "System.h"
25    
26     /**************************************************************/
27    
28     void Finley_SystemMatrix_nullifyRowsAndCols(Finley_SystemMatrix* A,
29     escriptDataC* mask_row, escriptDataC* mask_col, double main_diagonal_value) {
30    
31     /* check structure: */
32    
33     if (!isExpanded(mask_row)) {
34     Finley_ErrorCode=TYPE_ERROR;
35     sprintf(Finley_ErrorMsg,"row mask Data object has to be expanded");
36     } else if (!isExpanded(mask_col)) {
37     Finley_ErrorCode=TYPE_ERROR;
38     sprintf(Finley_ErrorMsg,"column mask Data object has to be expanded");
39     } else if (getLength(mask_col) != A ->col_block_size *A ->num_cols ) {
40     Finley_ErrorCode=TYPE_ERROR;
41     sprintf(Finley_ErrorMsg,"column mask vector and matrix column dimensions don't match.");
42     } else if (getLength(mask_row) != A ->row_block_size * A ->num_rows) {
43     Finley_ErrorCode=TYPE_ERROR;
44     sprintf(Finley_ErrorMsg,"row mask vector and matrix row dimensions don't match.");
45     }
46    
47     /* do the operation: */
48     if (Finley_ErrorCode==NO_ERROR) {
49     Finley_SystemMatrixNullify(A,getSampleData(mask_row,0),getSampleData(mask_col,0),main_diagonal_value);
50     }
51    
52     return;
53     }
54    
55     /****************************************************************/
56    
57     /* this is the actual ellimination */
58    
59     void Finley_SystemMatrixNullify(Finley_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
60    
61 jgs 102 maybelong ir,icol,iptr,icb,irb,irow,ic,l;
62 jgs 82
63     if (A ->col_block_size==1 && A ->row_block_size ==1) {
64     switch(A->type) {
65     case CSR:
66     #pragma omp parallel for private(irow, iptr,icol) schedule(static)
67 jgs 102 for (irow=0;irow< A->pattern->n_ptr;irow++) {
68 jgs 82 /* TODO: test mask_row here amd not inside every loop */
69 jgs 102 for (iptr=A->pattern->ptr[irow]-PTR_OFFSET;iptr<A->pattern->ptr[irow+1]-PTR_OFFSET; iptr++) {
70     icol=A->pattern->index[iptr]-INDEX_OFFSET;
71 jgs 82 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
72     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
73     A->val[iptr]=main_diagonal_value;
74     } else {
75     A->val[iptr]=0;
76     }
77     }
78     }
79     }
80     break;
81     case CSC:
82     #pragma omp parallel for private(irow, iptr,icol) schedule(static)
83 jgs 102 for (icol=0;icol< A->pattern->n_ptr;icol++) {
84     for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {
85     irow=A->pattern->index[iptr]-INDEX_OFFSET;
86 jgs 82 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
87     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
88     A->val[iptr]=main_diagonal_value;
89     } else {
90     A->val[iptr]=0;
91     }
92     }
93     }
94     }
95     break;
96     default:
97     Finley_ErrorCode=TYPE_ERROR;
98     sprintf(Finley_ErrorMsg,"Unknown matrix type.");
99     } /* switch A->type */
100     } else {
101     switch(A->type) {
102     case CSR:
103 jgs 102 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
104     for (ir=0;ir< A->pattern->n_ptr;ir++) {
105     for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
106 jgs 82 for (irb=0;irb< A->row_block_size;irb++) {
107     irow=irb+A->row_block_size*ir;
108     for (icb=0;icb< A->col_block_size;icb++) {
109 jgs 102 icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
110 jgs 82 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111 jgs 102 l=iptr*A->block_size+irb+A->row_block_size*icb;
112 jgs 82 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
113 jgs 102 A->val[l]=main_diagonal_value;
114 jgs 82 } else {
115 jgs 102 A->val[l]=0;
116 jgs 82 }
117     }
118     }
119     }
120     }
121     }
122     break;
123     case CSC:
124 jgs 102 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
125     for (ic=0;ic< A->pattern->n_ptr;ic++) {
126     for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
127 jgs 82 for (irb=0;irb< A->row_block_size;irb++) {
128 jgs 102 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
129 jgs 82 for (icb=0;icb< A->col_block_size;icb++) {
130     icol=icb+A->col_block_size*ic;
131     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
132 jgs 102 l=iptr*A->block_size+irb+A->row_block_size*icb;
133 jgs 82 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
134 jgs 102 A->val[l]=main_diagonal_value;
135 jgs 82 } else {
136 jgs 102 A->val[l]=0;
137 jgs 82 }
138     }
139     }
140     }
141     }
142     }
143     break;
144     default:
145     Finley_ErrorCode=TYPE_ERROR;
146     sprintf(Finley_ErrorMsg,"Unknown matrix type.");
147     } /* switch A->type */
148     }
149     return;
150     }
151    
152     /*
153     * $Log$
154 jgs 102 * Revision 1.4 2004/12/15 07:08:33 jgs
155 jgs 97 * *** empty log message ***
156 jgs 82 *
157 jgs 97 *
158 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26