/[escript]/trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c
ViewVC logotype

Annotation of /trunk/paso/src/SystemMatrix_nullifyRowsAndCols.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (hide annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 7 months ago) by ksteube
File MIME type: text/plain
File size: 4695 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


1 jgs 150 /* $Id$ */
2    
3 dhawcroft 631 /*
4     ********************************************************************************
5 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
6 dhawcroft 631 * *
7     * http://www.access.edu.au *
8     * Primary Business: Queensland, Australia *
9     * Licensed under the Open Software License version 3.0 *
10     * http://www.opensource.org/licenses/osl-3.0.php *
11     ********************************************************************************
12     */
13    
14 jgs 150 /**************************************************************/
15    
16     /* Paso: SystemMatrix */
17    
18     /* nullify rows and columns in the matrix */
19    
20     /* the rows and columns are marked by positive values in */
21     /* mask_row and mask_col. Values on the main diagonal */
22     /* which are marked to set to zero by both mask_row and */
23     /* mask_col are set to main_diagonal_value */
24    
25    
26     /**************************************************************/
27    
28     /* Copyrights by ACcESS Australia 2003 */
29     /* Author: gross@access.edu.au */
30    
31     /**************************************************************/
32    
33     #include "Paso.h"
34     #include "SystemMatrix.h"
35    
36     void Paso_SystemMatrix_nullifyRowsAndCols(Paso_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value) {
37    
38 gross 415 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
39 jgs 150 index_t ir,icol,iptr,icb,irb,irow,ic,l;
40    
41     if (A ->col_block_size==1 && A ->row_block_size ==1) {
42 gross 415 if (A->type & MATRIX_FORMAT_CSC) {
43 jgs 150 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
44 gross 415 for (icol=0;icol< A->pattern->n_ptr;icol++) {
45     for (iptr=A->pattern->ptr[icol]-index_offset;iptr<A->pattern->ptr[icol+1]-index_offset; iptr++) {
46     irow=A->pattern->index[iptr]-index_offset;
47 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
48     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
49     A->val[iptr]=main_diagonal_value;
50     } else {
51     A->val[iptr]=0;
52     }
53     }
54     }
55     }
56 gross 415 } else {
57 jgs 150 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
58 gross 415 for (irow=0;irow< A->pattern->n_ptr;irow++) {
59     /* TODO: test mask_row here amd not inside every loop */
60     for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
61     icol=A->pattern->index[iptr]-index_offset;
62 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
63     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
64     A->val[iptr]=main_diagonal_value;
65     } else {
66     A->val[iptr]=0;
67     }
68     }
69     }
70     }
71 gross 415 }
72 jgs 150 } else {
73 gross 415 if (A->type & MATRIX_FORMAT_CSC) {
74     #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
75     for (ic=0;ic< A->pattern->n_ptr;ic++) {
76     for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
77 jgs 150 for (irb=0;irb< A->row_block_size;irb++) {
78 gross 415 irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
79 jgs 150 for (icb=0;icb< A->col_block_size;icb++) {
80 gross 415 icol=icb+A->col_block_size*ic;
81 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
82     l=iptr*A->block_size+irb+A->row_block_size*icb;
83     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
84     A->val[l]=main_diagonal_value;
85     } else {
86     A->val[l]=0;
87     }
88     }
89     }
90     }
91     }
92     }
93 gross 415 } else {
94     #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
95     for (ir=0;ir< A->pattern->n_ptr;ir++) {
96     for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
97 jgs 150 for (irb=0;irb< A->row_block_size;irb++) {
98 gross 415 irow=irb+A->row_block_size*ir;
99 jgs 150 for (icb=0;icb< A->col_block_size;icb++) {
100 gross 415 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
101 jgs 150 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
102     l=iptr*A->block_size+irb+A->row_block_size*icb;
103     if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
104     A->val[l]=main_diagonal_value;
105     } else {
106     A->val[l]=0;
107     }
108     }
109     }
110     }
111     }
112     }
113 gross 415 }
114 jgs 150 }
115     return;
116     }
117    
118     /*
119     * $Log$
120     * Revision 1.2 2005/09/15 03:44:39 jgs
121     * Merge of development branch dev-02 back to main trunk on 2005-09-15
122     *
123     * Revision 1.1.2.1 2005/09/05 06:29:48 gross
124     * These files have been extracted from finley to define a stand alone libray for iterative
125     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
126     * has not been tested yet.
127     *
128     *
129     */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26