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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 6056 byte(s)
Don't panic.
Updating copyright stamps

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * 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
14
15 /**************************************************************/
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: Lutz Gross, l.gross@uq.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 #pragma ivdep
42 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 #pragma ivdep
61 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 index_t icol,iptr,icb,irb,irow,ic,l;
76 #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 #pragma ivdep
82 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 index_t ir,icol,iptr,icb,irb,irow,l;
100 #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 #pragma ivdep
106 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 #pragma ivdep
128 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 index_t ir,icol,iptr,icb,irb,irow,l;
142 #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 #pragma ivdep
149 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