/[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 100 - (hide annotations)
Wed Dec 15 03:48:48 2004 UTC (15 years, 5 months ago) by jgs
File MIME type: text/plain
File size: 5142 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 100 maybelong ir,icol,iptr,icb,irb,irow,ic;
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 100 for (irow=0;irow< A->num_rows;irow++) {
68 jgs 82 /* TODO: test mask_row here amd not inside every loop */
69 jgs 100 for (iptr=A->ptr[irow]-PTR_OFFSET;iptr<A->ptr[irow+1]-PTR_OFFSET; iptr++) {
70     icol=A->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 100 for (icol=0;icol< A->num_cols;icol++) {
84     for (iptr=A->ptr[icol]-PTR_OFFSET;iptr<A->ptr[icol+1]-PTR_OFFSET; iptr++) {
85     irow=A->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 100 #pragma omp parallel for private(irow, iptr,icol,ir,irb,icb) schedule(static)
104     for (ir=0;ir< A->num_rows;ir++) {
105     for (iptr=A->ptr[ir]-PTR_OFFSET;iptr<A->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 100 icol=icb+A->col_block_size*(A->index[iptr]-INDEX_OFFSET);
110 jgs 82 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
112 jgs 100 A->val[iptr]=main_diagonal_value;
113 jgs 82 } else {
114 jgs 100 A->val[iptr]=0;
115 jgs 82 }
116     }
117     }
118     }
119     }
120     }
121     break;
122     case CSC:
123 jgs 100 #pragma omp parallel for private(irow, iptr,icol,ic,irb,icb) schedule(static)
124     for (ic=0;ic< A->num_cols;ic++) {
125     for (iptr=A->ptr[ic]-PTR_OFFSET;iptr<A->ptr[ic+1]-PTR_OFFSET; iptr++) {
126 jgs 82 for (irb=0;irb< A->row_block_size;irb++) {
127 jgs 100 irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);
128 jgs 82 for (icb=0;icb< A->col_block_size;icb++) {
129     icol=icb+A->col_block_size*ic;
130     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
131     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
132 jgs 100 A->val[iptr]=main_diagonal_value;
133 jgs 82 } else {
134 jgs 100 A->val[iptr]=0;
135 jgs 82 }
136     }
137     }
138     }
139     }
140     }
141     break;
142     default:
143     Finley_ErrorCode=TYPE_ERROR;
144     sprintf(Finley_ErrorMsg,"Unknown matrix type.");
145     } /* switch A->type */
146     }
147     return;
148     }
149    
150     /*
151     * $Log$
152 jgs 100 * Revision 1.3 2004/12/15 03:48:46 jgs
153 jgs 97 * *** empty log message ***
154 jgs 82 *
155 jgs 97 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
156     * initial import of project esys2
157     *
158 jgs 82 * Revision 1.1 2004/07/02 04:21:13 gross
159     * Finley C code has been included
160     *
161     *
162     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26