/[escript]/trunk/paso/src/SparseMatrix_getSubmatrix.cpp
ViewVC logotype

Diff of /trunk/paso/src/SparseMatrix_getSubmatrix.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/paso/src/SparseMatrix_getSubmatrix.c revision 1315 by ksteube, Tue Sep 25 02:41:13 2007 UTC trunk/paso/src/SparseMatrix_getSubmatrix.cpp revision 4346 by jfenwick, Tue Apr 2 04:46:45 2013 UTC
# Line 1  Line 1 
1    
2  /* $Id: SparseMatrix_getSubmatrix.c 1306 2007-09-18 05:51:09Z ksteube $ */  /*****************************************************************************
3    *
4    * Copyright (c) 2003-2013 by University of Queensland
5    * http://www.uq.edu.au
6    *
7    * Primary Business: Queensland, Australia
8    * Licensed under the Open Software License version 3.0
9    * http://www.opensource.org/licenses/osl-3.0.php
10    *
11    * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
 /*******************************************************  
  *  
  *           Copyright 2003-2007 by ACceSS MNRF  
  *       Copyright 2007 by University of Queensland  
  *  
  *                http://esscc.uq.edu.au  
  *        Primary Business: Queensland, Australia  
  *  Licensed under the Open Software License version 3.0  
  *     http://www.opensource.org/licenses/osl-3.0.php  
  *  
  *******************************************************/  
16    
17  /**************************************************************/  /************************************************************************************/
18    
19  /* Paso: SparseMatrix */  /* Paso: SparseMatrix */
20    
21  /**************************************************************/  /************************************************************************************/
22    
23  /* Copyrights by ACcESS Australia 2003, 2004,2005 */  /* Copyrights by ACcESS Australia 2003, 2004,2005 */
24  /* Author: gross@access.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
25    
26  /**************************************************************/  /************************************************************************************/
27    
28  #include "Paso.h"  #include "Paso.h"
29  #include "SparseMatrix.h"  #include "SparseMatrix.h"
30  #include "PasoUtil.h"  #include "PasoUtil.h"
31    
32  /**************************************************************  /************************************************************************************
33    
34      returns the submatrix of A where rows are gathered by index row_list      Returns the submatrix of A where rows are gathered by index row_list
35      and columns are selected by non-negative values of new_col_index.      and columns are selected by non-negative values of new_col_index.
36      if new_col_index[i]>-1 new_col_index[i] gives the column of i in      If new_col_index[i]>-1 new_col_index[i] gives the column of i in
37      the returned submatrix      the returned submatrix.
   
38  */  */
39    
40    
# Line 44  Paso_SparseMatrix* Paso_SparseMatrix_get Line 44  Paso_SparseMatrix* Paso_SparseMatrix_get
44        index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);        index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
45        int i,k,tmp,m,subpattern_row;        int i,k,tmp,m,subpattern_row;
46        int type=A->type;        int type=A->type;
47        Paso_resetError();        Esys_resetError();
48        if (A->type & MATRIX_FORMAT_CSC) {        if (A->type & MATRIX_FORMAT_CSC) {
49            Paso_setError(TYPE_ERROR,"gathering submatrices supports CSR matrix format only.");            Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_getSubmatrix: gathering submatrices supports CSR matrix format only.");
50        } else {        } else {
51           sub_pattern=Paso_Pattern_getSubpattern(A->pattern,n_row_sub,n_col_sub,row_list,new_col_index);           sub_pattern=Paso_Pattern_getSubpattern(A->pattern,n_row_sub,n_col_sub,row_list,new_col_index);
52           if (Paso_noError()) {           if (Esys_noError()) {
53              /* create the return object */              /* create the return object */
54              out=Paso_SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size);              out=Paso_SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size,TRUE);
55              if (Paso_noError()) {              if (Esys_noError()) {
56                   #pragma omp parallel for private(i,k,m,subpattern_row,tmp) schedule(static)                   #pragma omp parallel for private(i,k,m,subpattern_row,tmp) schedule(static)
57                   for (i=0;i<n_row_sub;++i) {                   for (i=0;i<n_row_sub;++i) {
58                       subpattern_row=row_list[i];                       subpattern_row=row_list[i];
# Line 62  Paso_SparseMatrix* Paso_SparseMatrix_get Line 62  Paso_SparseMatrix* Paso_SparseMatrix_get
62                             #pragma ivdep                             #pragma ivdep
63                             for (m=out->pattern->ptr[i]-index_offset;m<out->pattern->ptr[i+1]-index_offset;++m) {                             for (m=out->pattern->ptr[i]-index_offset;m<out->pattern->ptr[i+1]-index_offset;++m) {
64                                 if (out->pattern->index[m]==tmp+index_offset) {                                 if (out->pattern->index[m]==tmp+index_offset) {
65                                     Paso_copyDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size]));                                     Paso_copyShortDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size]));
66                                     break;                                     break;
67                                 }                                 }
68                             }                             }
# Line 75  Paso_SparseMatrix* Paso_SparseMatrix_get Line 75  Paso_SparseMatrix* Paso_SparseMatrix_get
75        }        }
76        return out;        return out;
77  }  }
78    
79    Paso_SparseMatrix* Paso_SparseMatrix_getBlock(Paso_SparseMatrix* A, int blockid){
80          
81          dim_t blocksize=A->row_block_size,i;
82          dim_t n=A->numRows;
83          index_t iptr;
84          
85          Paso_SparseMatrix* out=NULL;
86        
87          out=Paso_SparseMatrix_alloc(A->type,A->pattern, 1, 1, 0);
88          
89          if (blocksize==1) {
90                if (blockid==1) {
91                      #pragma omp parallel for private(i,iptr) schedule(static)
92                      for(i=0;i<n;++i) {
93                          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
94                                  out->val[iptr]=A->val[iptr];
95                          }
96                      }
97                }
98                else {
99                      Esys_setError(VALUE_ERROR, "Paso_SparseMatrix_getBlock: Requested and actual block sizes do not match.");
100                }
101          }
102          else if (blocksize==2) {
103                if (blockid==1) {
104                      #pragma omp parallel for private(i,iptr) schedule(static)
105                      for(i=0;i<n;i++) {
106                            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
107                                  out->val[iptr]=A->val[4*iptr];
108                            }
109                      }
110                }
111                else if (blockid==2) {
112                      #pragma omp parallel for private(i,iptr) schedule(static)
113                      for(i=0;i<n;i++) {
114                            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
115                                  out->val[iptr]=A->val[4*iptr+3];
116                            }
117                      }
118                }
119                else {
120                      Esys_setError(VALUE_ERROR,"Paso_SparseMatrix_getBlock: Requested and actual block sizes do not match.");
121                }
122          }
123          else if (blocksize==3) {
124                if (blockid==1) {
125                      #pragma omp parallel for private(i,iptr) schedule(static)
126                      for(i=0;i<n;i++) {
127                            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
128                                        out->val[iptr]=A->val[9*iptr];
129                            }
130                      }
131                }
132                else if (blockid==2) {
133                      #pragma omp parallel for private(i,iptr) schedule(static)
134                      for(i=0;i<n;i++) {
135                            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
136                                        out->val[iptr]=A->val[9*iptr+4];
137                            }
138                      }
139                }
140                else if (blockid==3) {
141                      #pragma omp parallel for private(i,iptr) schedule(static)
142                      for(i=0;i<n;i++) {
143                            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
144                                        out->val[iptr]=A->val[9*iptr+8];
145                            }
146                      }
147                }
148                else {
149                      Esys_setError(VALUE_ERROR,"Paso_SparseMatrix_getBlock: Requested and actual block sizes do not match.");
150                }
151          }
152                
153             return out;
154    }

Legend:
Removed from v.1315  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26