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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 7 months ago) by ksteube
File MIME type: text/plain
File size: 5239 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


1 /* $Id$ */
2
3 /*
4 ********************************************************************************
5 * Copyright 2006 by ACcESS MNRF *
6 * *
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 /**************************************************************/
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 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
39 index_t ir,icol,iptr,icb,irb,irow,ic,l;
40
41 if (A ->col_block_size==1 && A ->row_block_size ==1) {
42 if (A->type & MATRIX_FORMAT_CSC) {
43 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
44 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 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 } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
57 /*
58 Plan: Collect non-zeros on each CPU,
59 use MPI to concatenate,
60 sort with qsort,
61 and then use binary search to check if a row is non-zero
62 */
63 fprintf(stderr, "SystemMatrix_nullifyRowsAndCols: need to implement MATRIX_FORMAT_TRILINOS_CRS (and for block sizes != 1)\n");
64 exit(1);
65 } else {
66 #pragma omp parallel for private(irow, iptr,icol) schedule(static)
67 for (irow=0;irow< A->pattern->n_ptr;irow++) {
68 /* TODO: test mask_row here amd not inside every loop */
69 for (iptr=A->pattern->ptr[irow]-index_offset;iptr<A->pattern->ptr[irow+1]-index_offset; iptr++) {
70 icol=A->pattern->index[iptr]-index_offset;
71 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
72 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
73 A->val[iptr]=main_diagonal_value;
74 } else {
75 A->val[iptr]=0;
76 }
77 }
78 }
79 }
80 }
81 } else {
82 if (A->type & MATRIX_FORMAT_CSC) {
83 #pragma omp parallel for private(l,irow, iptr,icol,ic,irb,icb) schedule(static)
84 for (ic=0;ic< A->pattern->n_ptr;ic++) {
85 for (iptr=A->pattern->ptr[ic]-index_offset;iptr<A->pattern->ptr[ic+1]-index_offset; iptr++) {
86 for (irb=0;irb< A->row_block_size;irb++) {
87 irow=irb+A->row_block_size*(A->pattern->index[iptr]-index_offset);
88 for (icb=0;icb< A->col_block_size;icb++) {
89 icol=icb+A->col_block_size*ic;
90 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
91 l=iptr*A->block_size+irb+A->row_block_size*icb;
92 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
93 A->val[l]=main_diagonal_value;
94 } else {
95 A->val[l]=0;
96 }
97 }
98 }
99 }
100 }
101 }
102 } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
103 fprintf(stderr, "SystemMatrix_nullifyRowsAndCols 2: need to implement MATRIX_FORMAT_TRILINOS_CRS\n");
104 exit(1);
105 } else {
106 #pragma omp parallel for private(l,irow, iptr,icol,ir,irb,icb) schedule(static)
107 for (ir=0;ir< A->pattern->n_ptr;ir++) {
108 for (iptr=A->pattern->ptr[ir]-index_offset;iptr<A->pattern->ptr[ir+1]-index_offset; iptr++) {
109 for (irb=0;irb< A->row_block_size;irb++) {
110 irow=irb+A->row_block_size*ir;
111 for (icb=0;icb< A->col_block_size;icb++) {
112 icol=icb+A->col_block_size*(A->pattern->index[iptr]-index_offset);
113 if (mask_col[icol]>0. || mask_row[irow]>0. ) {
114 l=iptr*A->block_size+irb+A->row_block_size*icb;
115 if (irow==icol && mask_col[icol]>0. && mask_row[irow]>0.) {
116 A->val[l]=main_diagonal_value;
117 } else {
118 A->val[l]=0;
119 }
120 }
121 }
122 }
123 }
124 }
125 }
126 }
127 return;
128 }
129
130 /*
131 * $Log$
132 * Revision 1.2 2005/09/15 03:44:39 jgs
133 * Merge of development branch dev-02 back to main trunk on 2005-09-15
134 *
135 * Revision 1.1.2.1 2005/09/05 06:29:48 gross
136 * These files have been extracted from finley to define a stand alone libray for iterative
137 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
138 * has not been tested yet.
139 *
140 *
141 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26