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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 7 months ago) by jgs
File MIME type: text/plain
File size: 5049 byte(s)
Initial revision

1 /* $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 maybelong ir,icol,iptr,icb,irb,irow,ic;
62
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 for (irow=0;irow< A->num_rows;irow++) {
68 /* TODO: test mask_row here amd not inside every loop */
69 for (iptr=A->ptr[irow]-PTR_OFFSET;iptr<A->ptr[irow+1]-PTR_OFFSET; iptr++) {
70 icol=A->index[iptr]-INDEX_OFFSET;
71 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 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 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 #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 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 icol=icb+A->col_block_size*(A->index[iptr]-INDEX_OFFSET);
110 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
112 A->val[iptr]=main_diagonal_value;
113 } else {
114 A->val[iptr]=0;
115 }
116 }
117 }
118 }
119 }
120 }
121 break;
122 case CSC:
123 #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 for (irb=0;irb< A->row_block_size;irb++) {
127 irow=irb+A->row_block_size*(A->index[iptr]-INDEX_OFFSET);
128 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 A->val[iptr]=main_diagonal_value;
133 } else {
134 A->val[iptr]=0;
135 }
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 * Revision 1.1 2004/10/26 06:53:57 jgs
153 * Initial revision
154 *
155 * Revision 1.1 2004/07/02 04:21:13 gross
156 * Finley C code has been included
157 *
158 *
159 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26