/[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 633 - (hide annotations)
Thu Mar 23 05:37:00 2006 UTC (13 years, 5 months ago) by dhawcroft
File MIME type: text/plain
File size: 4695 byte(s)


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