/[escript]/trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c
ViewVC logotype

Contents of /trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 10 months ago) by jgs
File MIME type: text/plain
File size: 4258 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Paso: 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 "Paso.h"
23 #include "SystemMatrix.h"
24
25 void Paso_SystemMatrix_nullifyRowsAndCols(Paso_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
26
27 index_t ir,icol,iptr,icb,irb,irow,ic,l;
28
29 if (A ->col_block_size==1 && A ->row_block_size ==1) {
30 switch(A->type) {
31 case CSR:
32 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
33 for (irow=0;irow< A->pattern->n_ptr;irow++) {
34 /* TODO: test mask_row here amd not inside every loop */
35 for (iptr=A->pattern->ptr[irow]-PTR_OFFSET;iptr<A->pattern->ptr[irow+1]-PTR_OFFSET; iptr++) {
36 icol=A->pattern->index[iptr]-INDEX_OFFSET;
37 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
38 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
39 A->val[iptr]=main_diagonal_value;
40 } else {
41 A->val[iptr]=0;
42 }
43 }
44 }
45 }
46 break;
47 case CSC:
48 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
49 for (icol=0;icol< A->pattern->n_ptr;icol++) {
50 for (iptr=A->pattern->ptr[icol]-PTR_OFFSET;iptr<A->pattern->ptr[icol+1]-PTR_OFFSET; iptr++) {
51 irow=A->pattern->index[iptr]-INDEX_OFFSET;
52 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
53 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
54 A->val[iptr]=main_diagonal_value;
55 } else {
56 A->val[iptr]=0;
57 }
58 }
59 }
60 }
61 break;
62 default:
63 Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");
64 } /* switch A->type */
65 } else {
66 switch(A->type) {
67 case CSR:
68 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
69 for (ir=0;ir< A->pattern->n_ptr;ir++) {
70 for (iptr=A->pattern->ptr[ir]-PTR_OFFSET;iptr<A->pattern->ptr[ir+1]-PTR_OFFSET; iptr++) {
71 for (irb=0;irb< A->row_block_size;irb++) {
72 irow=irb+A->row_block_size*ir;
73 for (icb=0;icb< A->col_block_size;icb++) {
74 icol=icb+A->col_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
75 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
76 l=iptr*A->block_size+irb+A->row_block_size*icb;
77 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
78 A->val[l]=main_diagonal_value;
79 } else {
80 A->val[l]=0;
81 }
82 }
83 }
84 }
85 }
86 }
87 break;
88 case CSC:
89 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
90 for (ic=0;ic< A->pattern->n_ptr;ic++) {
91 for (iptr=A->pattern->ptr[ic]-PTR_OFFSET;iptr<A->pattern->ptr[ic+1]-PTR_OFFSET; iptr++) {
92 for (irb=0;irb< A->row_block_size;irb++) {
93 irow=irb+A->row_block_size*(A->pattern->index[iptr]-INDEX_OFFSET);
94 for (icb=0;icb< A->col_block_size;icb++) {
95 icol=icb+A->col_block_size*ic;
96 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
97 l=iptr*A->block_size+irb+A->row_block_size*icb;
98 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
99 A->val[l]=main_diagonal_value;
100 } else {
101 A->val[l]=0;
102 }
103 }
104 }
105 }
106 }
107 }
108 break;
109 default:
110 Paso_setError(TYPE_ERROR,"Unknown matrix type in row and column elimination.");
111 } /* switch A->type */
112 }
113 return;
114 }
115
116 /*
117 * $Log$
118 * Revision 1.2 2005/09/15 03:44:39 jgs
119 * Merge of development branch dev-02 back to main trunk on 2005-09-15
120 *
121 * Revision 1.1.2.1 2005/09/05 06:29:48 gross
122 * These files have been extracted from finley to define a stand alone libray for iterative
123 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
124 * has not been tested yet.
125 *
126 *
127 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26