/[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

revision 4828 by caltinay, Tue Apr 1 03:50:23 2014 UTC revision 4829 by caltinay, Thu Apr 3 04:02:53 2014 UTC
# Line 41  namespace paso { Line 41  namespace paso {
41  */  */
42    
43    
44  SparseMatrix* SparseMatrix_getSubmatrix(const SparseMatrix* A,int n_row_sub,  SparseMatrix_ptr SparseMatrix::getSubmatrix(int n_row_sub, int n_col_sub,
45                                          int n_col_sub, const index_t* row_list,                                              const index_t* row_list,
46                                          const index_t* new_col_index)                                              const index_t* new_col_index) const
47  {  {
48      SparseMatrix* out=NULL;      SparseMatrix_ptr out;
     index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);  
     int i,k,tmp,m,subpattern_row;  
     int type=A->type;  
49      Esys_resetError();      Esys_resetError();
50      if (A->type & MATRIX_FORMAT_CSC) {      if (type & MATRIX_FORMAT_CSC) {
51          Esys_setError(TYPE_ERROR, "SparseMatrix_getSubmatrix: gathering submatrices supports CSR matrix format only.");          Esys_setError(TYPE_ERROR, "SparseMatrix::getSubmatrix: gathering submatrices supports CSR matrix format only.");
52      } else {          return out;
53          Pattern_ptr sub_pattern(A->pattern->getSubpattern(n_row_sub, n_col_sub,      }
54                                            row_list, new_col_index));  
55        const index_t index_offset = (type & MATRIX_FORMAT_OFFSET1 ? 1:0);
56        Pattern_ptr sub_pattern(pattern->getSubpattern(n_row_sub, n_col_sub,
57                                                       row_list, new_col_index));
58        if (Esys_noError()) {
59            // create the return object
60            out.reset(new SparseMatrix(type, sub_pattern, row_block_size,
61                                       col_block_size, true));
62          if (Esys_noError()) {          if (Esys_noError()) {
63              /* create the return object */  #pragma omp parallel for
64              out=SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size,TRUE);              for (int i=0; i<n_row_sub; ++i) {
65              if (Esys_noError()) {                  const index_t subpattern_row = row_list[i];
66                  #pragma omp parallel for private(i,k,m,subpattern_row,tmp) schedule(static)                  for (int k=pattern->ptr[subpattern_row]-index_offset;
67                  for (i=0;i<n_row_sub;++i) {                          k < pattern->ptr[subpattern_row+1]-index_offset; ++k) {
68                      subpattern_row=row_list[i];                      index_t tmp=new_col_index[pattern->index[k]-index_offset];
69                      for (k=A->pattern->ptr[subpattern_row]-index_offset;k<A->pattern->ptr[subpattern_row+1]-index_offset;++k) {                      if (tmp > -1) {
70                          tmp=new_col_index[A->pattern->index[k]-index_offset];                          #pragma ivdep
71                          if (tmp>-1) {                          for (index_t m=out->pattern->ptr[i]-index_offset;
72                              #pragma ivdep                                  m < out->pattern->ptr[i+1]-index_offset; ++m) {
73                              for (m=out->pattern->ptr[i]-index_offset;m<out->pattern->ptr[i+1]-index_offset;++m) {                              if (out->pattern->index[m]==tmp+index_offset) {
74                                  if (out->pattern->index[m]==tmp+index_offset) {                                  Paso_copyShortDouble(block_size, &val[k*block_size], &out->val[m*block_size]);
75                                     Paso_copyShortDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size]));                                  break;
                                    break;  
                                 }  
76                              }                              }
77                          }                          }
78                      }                      }
# Line 81  SparseMatrix* SparseMatrix_getSubmatrix( Line 83  SparseMatrix* SparseMatrix_getSubmatrix(
83      return out;      return out;
84  }  }
85    
86  SparseMatrix* SparseMatrix_getBlock(const SparseMatrix* A, int blockid)  SparseMatrix_ptr SparseMatrix::getBlock(int blockid) const
87  {  {
88      dim_t blocksize=A->row_block_size,i;      const dim_t blocksize = row_block_size;
89      dim_t n=A->numRows;      const dim_t n = numRows;
90      index_t iptr;      SparseMatrix_ptr out(new SparseMatrix(type, pattern, 1, 1, 0));
     SparseMatrix* out = SparseMatrix_alloc(A->type, A->pattern, 1, 1, 0);  
91    
92      if (blocksize==1) {      if (blocksize==1) {
93          if (blockid==1) {          if (blockid==1) {
94              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
95              for(i=0;i<n;++i) {              for (dim_t i=0; i<n; ++i) {
96                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
97                      out->val[iptr]=A->val[iptr];                      out->val[iptr] = val[iptr];
98                  }                  }
99              }              }
100          } else {          } else {
101              Esys_setError(VALUE_ERROR, "SparseMatrix_getBlock: Requested and actual block sizes do not match.");              Esys_setError(VALUE_ERROR, "SparseMatrix::getBlock: Invalid block ID requested.");
102          }          }
103      } else if (blocksize==2) {      } else if (blocksize==2) {
104          if (blockid==1) {          if (blockid==1) {
105              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
106              for(i=0;i<n;i++) {              for (dim_t i=0; i<n; i++) {
107                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
108                      out->val[iptr]=A->val[4*iptr];                      out->val[iptr] = val[4*iptr];
109                  }                  }
110              }              }
111          } else if (blockid==2) {          } else if (blockid==2) {
112              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
113              for(i=0;i<n;i++) {              for (dim_t i=0; i<n; i++) {
114                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
115                      out->val[iptr]=A->val[4*iptr+3];                      out->val[iptr] = val[4*iptr+3];
116                  }                  }
117              }              }
118          } else {          } else {
119              Esys_setError(VALUE_ERROR,"SparseMatrix_getBlock: Requested and actual block sizes do not match.");              Esys_setError(VALUE_ERROR, "SparseMatrix::getBlock: Invalid block ID requested.");
120          }          }
121      } else if (blocksize==3) {      } else if (blocksize==3) {
122          if (blockid==1) {          if (blockid==1) {
123              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
124              for(i=0;i<n;i++) {              for (dim_t i=0; i<n; i++) {
125                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
126                      out->val[iptr]=A->val[9*iptr];                      out->val[iptr] = val[9*iptr];
127                  }                  }
128              }              }
129          } else if (blockid==2) {          } else if (blockid==2) {
130              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
131              for(i=0;i<n;i++) {              for (dim_t i=0; i<n; i++) {
132                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
133                      out->val[iptr]=A->val[9*iptr+4];                      out->val[iptr] = val[9*iptr+4];
134                  }                  }
135              }              }
136          } else if (blockid==3) {          } else if (blockid==3) {
137              #pragma omp parallel for private(i,iptr) schedule(static)  #pragma omp parallel for
138              for(i=0;i<n;i++) {              for (dim_t i=0; i<n; i++) {
139                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {                  for (index_t iptr=pattern->ptr[i]; iptr<pattern->ptr[i+1]; ++iptr) {
140                      out->val[iptr]=A->val[9*iptr+8];                      out->val[iptr] = val[9*iptr+8];
141                  }                  }
142              }              }
143          } else {          } else {
144              Esys_setError(VALUE_ERROR,"SparseMatrix_getBlock: Requested and actual block sizes do not match.");              Esys_setError(VALUE_ERROR, "SparseMatrix::getBlock: Invalid block ID requested.");
145          }          }
146      }      }
   
147      return out;      return out;
148  }  }
149    

Legend:
Removed from v.4828  
changed lines
  Added in v.4829

  ViewVC Help
Powered by ViewVC 1.1.26