/[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 1315 - (show annotations)
Tue Sep 25 02:41:13 2007 UTC (11 years, 11 months ago) by ksteube
File MIME type: text/plain
File size: 6065 byte(s)
Copied more files from MPI branch to trunk

1
2 /* $Id: SparseMatrix_nullifyRowsAndCols.c 1306 2007-09-18 05:51:09Z ksteube $ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* Paso: SparseMatrix */
19
20 /* nullify rows and columns in the matrix */
21
22 /* the rows and columns are marked by positive values in */
23 /* mask_row and mask_col. Values on the main diagonal */
24 /* which are marked to set to zero by both mask_row and */
25 /* mask_col are set to main_diagonal_value */
26
27
28 /**************************************************************/
29
30 /* Copyrights by ACcESS Australia 2003, 2007 */
31 /* Author: gross@access.edu.au */
32
33 /**************************************************************/
34
35 #include "Paso.h"
36 #include "SparseMatrix.h"
37
38 void Paso_SparseMatrix_nullifyRowsAndCols_CSC_BLK1(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
39 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
40 index_t irow, iptr, icol;
41 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
42 for (icol=0;icol< A->pattern->numOutput;icol++) {
43 for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
44 irow=A->pattern->index[iptr]-index_offset;
45 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
46 if (irow==icol) {
47 A->val[iptr]=main_diagonal_value;
48 } else {
49 A->val[iptr]=0;
50 }
51 }
52 }
53 }
54 }
55 void Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
56 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
57 index_t irow, iptr, icol;
58 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
59 for (irow=0;irow< A->pattern->numOutput;irow++) {
60 /* TODO: test mask_row here amd not inside every loop */
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 ir,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 for (icb=0;icb< A->col_block_size;icb++) {
82 icol=icb+A->col_block_size*ic;
83 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
84 l=iptr*A->block_size+irb+A->row_block_size*icb;
85 if (irow==icol) {
86 A->val[l]=main_diagonal_value;
87 } else {
88 A->val[l]=0;
89 }
90 }
91 }
92 }
93 }
94 }
95 }
96 void Paso_SparseMatrix_nullifyRowsAndCols_CSR(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
97 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
98 index_t ir,icol,iptr,icb,irb,irow,ic,l;
99 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
100 for (ir=0;ir< A->pattern->numOutput;ir++) {
101 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
102 for (irb=0;irb< A->row_block_size;irb++) {
103 irow=irb+A->row_block_size*ir;
104 for (icb=0;icb< A->col_block_size;icb++) {
105 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
106 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
107 l=iptr*A->block_size+irb+A->row_block_size*icb;
108 if (irow==icol) {
109 A->val[l]=main_diagonal_value;
110 } else {
111 A->val[l]=0;
112 }
113 }
114 }
115 }
116 }
117 }
118 }
119 void Paso_SparseMatrix_nullifyRows_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
120 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
121 index_t irow, iptr, icol;
122 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
123 for (irow=0;irow< A->pattern->numOutput;irow++) {
124 if (mask_row[irow]>0.) {
125 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
126 icol=A->pattern->index[iptr]-index_offset;
127 if (irow==icol) {
128 A->val[iptr]=main_diagonal_value;
129 } else {
130 A->val[iptr]=0;
131 }
132 }
133 }
134 }
135 }
136 void Paso_SparseMatrix_nullifyRows_CSR(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
137 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
138 index_t ir,icol,iptr,icb,irb,irow,ic,l;
139 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
140 for (ir=0;ir< A->pattern->numOutput;ir++) {
141 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
142 for (irb=0;irb< A->row_block_size;irb++) {
143 irow=irb+A->row_block_size*ir;
144 if (mask_row[irow]>0. ) {
145 for (icb=0;icb< A->col_block_size;icb++) {
146 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
147 l=iptr*A->block_size+irb+A->row_block_size*icb;
148 if (irow==icol) {
149 A->val[l]=main_diagonal_value;
150 } else {
151 A->val[l]=0;
152 }
153 }
154 }
155 }
156 }
157 }
158 }

  ViewVC Help
Powered by ViewVC 1.1.26