/[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 1563 - (show annotations)
Wed May 21 15:27:34 2008 UTC (12 years, 5 months ago) by gross
File MIME type: text/plain
File size: 6248 byte(s)
Some stuff added to make the FCT solver work with Dirichlet Boundary conditions for MPI. 
and a memory leak in the coupler fixed.


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 #pragma ivdep
44 for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
45 irow=A->pattern->index[iptr]-index_offset;
46 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
47 if (irow==icol) {
48 A->val[iptr]=main_diagonal_value;
49 } else {
50 A->val[iptr]=0;
51 }
52 }
53 }
54 }
55 }
56 void Paso_SparseMatrix_nullifyRowsAndCols_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
57 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
58 index_t irow, iptr, icol;
59 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
60 for (irow=0;irow< A->pattern->numOutput;irow++) {
61 /* TODO: test mask_row here amd not inside every loop */
62 #pragma ivdep
63 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
64 icol=A->pattern->index[iptr]-index_offset;
65 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
66 printf("eliminated %d %d\n",icol,irow);
67 if (irow==icol) {
68 A->val[iptr]=main_diagonal_value;
69 } else {
70 A->val[iptr]=0;
71 }
72 }
73 }
74 }
75 }
76 void Paso_SparseMatrix_nullifyRowsAndCols_CSC(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
77 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
78 index_t ir,icol,iptr,icb,irb,irow,ic,l;
79 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
80 for (ic=0;ic< A->pattern->numOutput;ic++) {
81 for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
82 for (irb=0;irb< A->row_block_size;irb++) {
83 irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
84 #pragma ivdep
85 for (icb=0;icb< A->col_block_size;icb++) {
86 icol=icb+A->col_block_size*ic;
87 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
88 l=iptr*A->block_size+irb+A->row_block_size*icb;
89 if (irow==icol) {
90 A->val[l]=main_diagonal_value;
91 } else {
92 A->val[l]=0;
93 }
94 }
95 }
96 }
97 }
98 }
99 }
100 void Paso_SparseMatrix_nullifyRowsAndCols_CSR(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
101 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
102 index_t ir,icol,iptr,icb,irb,irow,ic,l;
103 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
104 for (ir=0;ir< A->pattern->numOutput;ir++) {
105 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
106 for (irb=0;irb< A->row_block_size;irb++) {
107 irow=irb+A->row_block_size*ir;
108 #pragma ivdep
109 for (icb=0;icb< A->col_block_size;icb++) {
110 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
111 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
112 l=iptr*A->block_size+irb+A->row_block_size*icb;
113 if (irow==icol) {
114 A->val[l]=main_diagonal_value;
115 } else {
116 A->val[l]=0;
117 }
118 }
119 }
120 }
121 }
122 }
123 }
124 void Paso_SparseMatrix_nullifyRows_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
125 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
126 index_t irow, iptr, icol;
127 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
128 for (irow=0;irow< A->pattern->numOutput;irow++) {
129 if (mask_row[irow]>0.) {
130 #pragma ivdep
131 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
132 icol=A->pattern->index[iptr]-index_offset;
133 if (irow==icol) {
134 A->val[iptr]=main_diagonal_value;
135 } else {
136 A->val[iptr]=0;
137 }
138 }
139 }
140 }
141 }
142 void Paso_SparseMatrix_nullifyRows_CSR(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
143 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
144 index_t ir,icol,iptr,icb,irb,irow,ic,l;
145 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
146 for (ir=0;ir< A->pattern->numOutput;ir++) {
147 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
148 for (irb=0;irb< A->row_block_size;irb++) {
149 irow=irb+A->row_block_size*ir;
150 if (mask_row[irow]>0. ) {
151 #pragma ivdep
152 for (icb=0;icb< A->col_block_size;icb++) {
153 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
154 l=iptr*A->block_size+irb+A->row_block_size*icb;
155 if (irow==icol) {
156 A->val[l]=main_diagonal_value;
157 } else {
158 A->val[l]=0;
159 }
160 }
161 }
162 }
163 }
164 }
165 }

  ViewVC Help
Powered by ViewVC 1.1.26