/[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 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 6 months ago) by jgs
File MIME type: text/plain
File size: 5230 byte(s)
*** empty log message ***

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,l;
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->pattern->n_ptr;irow++) {
68 /* TODO: test mask_row here amd not inside every loop */
69 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 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->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 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(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 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->pattern->index[iptr]-INDEX_OFFSET);
110 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111 l=iptr*A->block_size+irb+A->row_block_size*icb;
112 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
113 A->val[l]=main_diagonal_value;
114 } else {
115 A->val[l]=0;
116 }
117 }
118 }
119 }
120 }
121 }
122 break;
123 case CSC:
124 #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 for (irb=0;irb< A->row_block_size;irb++) {
128 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
129 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 l=iptr*A->block_size+irb+A->row_block_size*icb;
133 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
134 A->val[l]=main_diagonal_value;
135 } else {
136 A->val[l]=0;
137 }
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 * Revision 1.4 2004/12/15 07:08:33 jgs
155 * *** empty log message ***
156 *
157 *
158 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26