/[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 1564 - (show annotations)
Thu May 22 09:31:33 2008 UTC (11 years, 5 months ago) by gross
File MIME type: text/plain
File size: 6208 byte(s)
some openmp dynamic scheduling for MVM.
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 if (irow==icol) {
67 A->val[iptr]=main_diagonal_value;
68 } else {
69 A->val[iptr]=0;
70 }
71 }
72 }
73 }
74 }
75 void Paso_SparseMatrix_nullifyRowsAndCols_CSC(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
76 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
77 index_t ir,icol,iptr,icb,irb,irow,ic,l;
78 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
79 for (ic=0;ic< A->pattern->numOutput;ic++) {
80 for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
81 for (irb=0;irb< A->row_block_size;irb++) {
82 irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
83 #pragma ivdep
84 for (icb=0;icb< A->col_block_size;icb++) {
85 icol=icb+A->col_block_size*ic;
86 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
87 l=iptr*A->block_size+irb+A->row_block_size*icb;
88 if (irow==icol) {
89 A->val[l]=main_diagonal_value;
90 } else {
91 A->val[l]=0;
92 }
93 }
94 }
95 }
96 }
97 }
98 }
99 void Paso_SparseMatrix_nullifyRowsAndCols_CSR(Paso_SparseMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
100 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
101 index_t ir,icol,iptr,icb,irb,irow,ic,l;
102 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
103 for (ir=0;ir< A->pattern->numOutput;ir++) {
104 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
105 for (irb=0;irb< A->row_block_size;irb++) {
106 irow=irb+A->row_block_size*ir;
107 #pragma ivdep
108 for (icb=0;icb< A->col_block_size;icb++) {
109 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
110 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
111 l=iptr*A->block_size+irb+A->row_block_size*icb;
112 if (irow==icol) {
113 A->val[l]=main_diagonal_value;
114 } else {
115 A->val[l]=0;
116 }
117 }
118 }
119 }
120 }
121 }
122 }
123 void Paso_SparseMatrix_nullifyRows_CSR_BLK1(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
124 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
125 index_t irow, iptr, icol;
126 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
127 for (irow=0;irow< A->pattern->numOutput;irow++) {
128 if (mask_row[irow]>0.) {
129 #pragma ivdep
130 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
131 icol=A->pattern->index[iptr]-index_offset;
132 if (irow==icol) {
133 A->val[iptr]=main_diagonal_value;
134 } else {
135 A->val[iptr]=0;
136 }
137 }
138 }
139 }
140 }
141 void Paso_SparseMatrix_nullifyRows_CSR(Paso_SparseMatrix* A, double* mask_row, double main_diagonal_value) {
142 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
143 index_t ir,icol,iptr,icb,irb,irow,ic,l;
144 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
145 for (ir=0;ir< A->pattern->numOutput;ir++) {
146 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
147 for (irb=0;irb< A->row_block_size;irb++) {
148 irow=irb+A->row_block_size*ir;
149 if (mask_row[irow]>0. ) {
150 #pragma ivdep
151 for (icb=0;icb< A->col_block_size;icb++) {
152 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
153 l=iptr*A->block_size+irb+A->row_block_size*icb;
154 if (irow==icol) {
155 A->val[l]=main_diagonal_value;
156 } else {
157 A->val[l]=0;
158 }
159 }
160 }
161 }
162 }
163 }
164 }

  ViewVC Help
Powered by ViewVC 1.1.26