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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 6046 byte(s)
Updating copyright notices
1 ksteube 1315
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1315
14 ksteube 1811
15 ksteube 1315 /**************************************************************/
16    
17     /* Paso: SparseMatrix */
18    
19     /* nullify rows and columns in the matrix */
20    
21     /* the rows and columns are marked by positive values in */
22     /* mask_row and mask_col. Values on the main diagonal */
23     /* which are marked to set to zero by both mask_row and */
24     /* mask_col are set to main_diagonal_value */
25    
26    
27     /**************************************************************/
28    
29     /* Author: gross@access.edu.au */
30    
31     /**************************************************************/
32    
33     #include "Paso.h"
34     #include "SparseMatrix.h"
35    
36     void Paso_SparseMatrix_nullifyRowsAndCols_CSC_BLK1(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
37     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38     index_t irow, iptr, icol;
39     #pragma omp parallel for private(irow, iptr,icol) schedule(static)
40     for (icol=0;icol< A->pattern->numOutput;icol++) {
41 gross 1552 #pragma ivdep
42 ksteube 1315 for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
43     irow=A->pattern->index[iptr]-index_offset;
44     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
45     if (irow==icol) {
46     A->val[iptr]=main_diagonal_value;
47     } else {
48     A->val[iptr]=0;
49     }
50     }
51     }
52     }
53     }
54     void Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
55     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
56     index_t irow, iptr, icol;
57     #pragma omp parallel for private(irow, iptr,icol) schedule(static)
58     for (irow=0;irow< A->pattern->numOutput;irow++) {
59     /* TODO: test mask_row here amd not inside every loop */
60 gross 1552 #pragma ivdep
61 ksteube 1315 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
62     icol=A->pattern->index[iptr]-index_offset;
63     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
64     if (irow==icol) {
65     A->val[iptr]=main_diagonal_value;
66     } else {
67     A->val[iptr]=0;
68     }
69     }
70     }
71     }
72     }
73     void Paso_SparseMatrix_nullifyRowsAndCols_CSC(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
74     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
75 phornby 1628 index_t icol,iptr,icb,irb,irow,ic,l;
76 ksteube 1315 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
77     for (ic=0;ic< A->pattern->numOutput;ic++) {
78     for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
79     for (irb=0;irb< A->row_block_size;irb++) {
80     irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
81 gross 1552 #pragma ivdep
82 ksteube 1315 for (icb=0;icb< A->col_block_size;icb++) {
83     icol=icb+A->col_block_size*ic;
84     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
85     l=iptr*A->block_size+irb+A->row_block_size*icb;
86     if (irow==icol) {
87     A->val[l]=main_diagonal_value;
88     } else {
89     A->val[l]=0;
90     }
91     }
92     }
93     }
94     }
95     }
96     }
97     void Paso_SparseMatrix_nullifyRowsAndCols_CSR(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
98     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
99 phornby 1628 index_t ir,icol,iptr,icb,irb,irow,l;
100 ksteube 1315 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
101     for (ir=0;ir< A->pattern->numOutput;ir++) {
102     for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
103     for (irb=0;irb< A->row_block_size;irb++) {
104     irow=irb+A->row_block_size*ir;
105 gross 1552 #pragma ivdep
106 ksteube 1315 for (icb=0;icb< A->col_block_size;icb++) {
107     icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
108     if (mask_col[icol]>0. || mask_row[irow]>0. ) {
109     l=iptr*A->block_size+irb+A->row_block_size*icb;
110     if (irow==icol) {
111     A->val[l]=main_diagonal_value;
112     } else {
113     A->val[l]=0;
114     }
115     }
116     }
117     }
118     }
119     }
120     }
121     void Paso_SparseMatrix_nullifyRows_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
122     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
123     index_t irow, iptr, icol;
124     #pragma omp parallel for private(irow, iptr,icol) schedule(static)
125     for (irow=0;irow< A->pattern->numOutput;irow++) {
126     if (mask_row[irow]>0.) {
127 gross 1552 #pragma ivdep
128 ksteube 1315 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
129     icol=A->pattern->index[iptr]-index_offset;
130     if (irow==icol) {
131     A->val[iptr]=main_diagonal_value;
132     } else {
133     A->val[iptr]=0;
134     }
135     }
136     }
137     }
138     }
139     void Paso_SparseMatrix_nullifyRows_CSR(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
140     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
141 phornby 1628 index_t ir,icol,iptr,icb,irb,irow,l;
142 ksteube 1315 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
143     for (ir=0;ir< A->pattern->numOutput;ir++) {
144     for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
145     for (irb=0;irb< A->row_block_size;irb++) {
146     irow=irb+A->row_block_size*ir;
147     if (mask_row[irow]>0. ) {
148 gross 1552 #pragma ivdep
149 ksteube 1315 for (icb=0;icb< A->col_block_size;icb++) {
150     icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
151     l=iptr*A->block_size+irb+A->row_block_size*icb;
152     if (irow==icol) {
153     A->val[l]=main_diagonal_value;
154     } else {
155     A->val[l]=0;
156     }
157     }
158     }
159     }
160     }
161     }
162     }

  ViewVC Help
Powered by ViewVC 1.1.26