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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/paso/src/SystemMatrix_nullifyRowsAndCols.c
File MIME type: text/plain
File size: 4258 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 150 /* $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