/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (12 years, 1 month ago) by phornby
Original Path: temp_trunk_copy/paso/src/SparseMatrix_nullifyRowsAndCols.c
File MIME type: text/plain
File size: 6065 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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     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