/[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 415 - (hide annotations)
Wed Jan 4 05:37:33 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/plain
File size: 4050 byte(s)
a better way of representing the matrix format type is introduced. this is needed for the Paradiso and UMFPACK interface
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 gross 415 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
28 jgs 150 index_t ir,icol,iptr,icb,irb,irow,ic,l;
29    
30     if (A ->col_block_size==1 && A ->row_block_size ==1) {
31 gross 415 if (A->type & MATRIX_FORMAT_CSC) {
32 jgs 150 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
33 gross 415 for (icol=0;icol< A->pattern->n_ptr;icol++) {
34     for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
35     irow=A->pattern->index[iptr]-index_offset;
36 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
37     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
38     A->val[iptr]=main_diagonal_value;
39     } else {
40     A->val[iptr]=0;
41     }
42     }
43     }
44     }
45 gross 415 } else {
46 jgs 150 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
47 gross 415 for (irow=0;irow< A->pattern->n_ptr;irow++) {
48     /* TODO: test mask_row here amd not inside every loop */
49     for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
50     icol=A->pattern->index[iptr]-index_offset;
51 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
52     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
53     A->val[iptr]=main_diagonal_value;
54     } else {
55     A->val[iptr]=0;
56     }
57     }
58     }
59     }
60 gross 415 }
61 jgs 150 } else {
62 gross 415 if (A->type & MATRIX_FORMAT_CSC) {
63     #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
64     for (ic=0;ic< A->pattern->n_ptr;ic++) {
65     for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
66 jgs 150 for (irb=0;irb< A->row_block_size;irb++) {
67 gross 415 irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
68 jgs 150 for (icb=0;icb< A->col_block_size;icb++) {
69 gross 415 icol=icb+A->col_block_size*ic;
70 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
71     l=iptr*A->block_size+irb+A->row_block_size*icb;
72     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
73     A->val[l]=main_diagonal_value;
74     } else {
75     A->val[l]=0;
76     }
77     }
78     }
79     }
80     }
81     }
82 gross 415 } else {
83     #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
84     for (ir=0;ir< A->pattern->n_ptr;ir++) {
85     for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
86 jgs 150 for (irb=0;irb< A->row_block_size;irb++) {
87 gross 415 irow=irb+A->row_block_size*ir;
88 jgs 150 for (icb=0;icb< A->col_block_size;icb++) {
89 gross 415 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
90 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
91     l=iptr*A->block_size+irb+A->row_block_size*icb;
92     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
93     A->val[l]=main_diagonal_value;
94     } else {
95     A->val[l]=0;
96     }
97     }
98     }
99     }
100     }
101     }
102 gross 415 }
103 jgs 150 }
104     return;
105     }
106    
107     /*
108     * $Log$
109     * Revision 1.2 2005/09/15 03:44:39 jgs
110     * Merge of development branch dev-02 back to main trunk on 2005-09-15
111     *
112     * Revision 1.1.2.1 2005/09/05 06:29:48 gross
113     * These files have been extracted from finley to define a stand alone libray for iterative
114     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
115     * has not been tested yet.
116     *
117     *
118     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26