/[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 1563 - (hide 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 ksteube 1315
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 gross 1552 #pragma ivdep
44 ksteube 1315 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 gross 1552 #pragma ivdep
63 ksteube 1315 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 gross 1563 printf("eliminated %d %d\n",icol,irow);
67 ksteube 1315 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 gross 1552 #pragma ivdep
85 ksteube 1315 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 gross 1552 #pragma ivdep
109 ksteube 1315 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 gross 1552 #pragma ivdep
131 ksteube 1315 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 gross 1552 #pragma ivdep
152 ksteube 1315 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