# Contents of /trunk/paso/src/SparseMatrix_getSubmatrix.cpp

Revision 4257 - (show annotations)
Wed Feb 27 03:42:40 2013 UTC (6 years, 7 months ago) by jfenwick
Original Path: branches/doubleplusgood/paso/src/SparseMatrix_getSubmatrix.c
File MIME type: text/plain
File size: 6308 byte(s)
```Some simple experiments for c++ Finley

```
 1 2 /***************************************************************************** 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 16 17 /************************************************************************************/ 18 19 /* Paso: SparseMatrix */ 20 21 /************************************************************************************/ 22 23 /* Copyrights by ACcESS Australia 2003, 2004,2005 */ 24 /* Author: Lutz Gross, l.gross@uq.edu.au */ 25 26 /************************************************************************************/ 27 28 #include "Paso.h" 29 #include "SparseMatrix.h" 30 #include "PasoUtil.h" 31 32 /************************************************************************************ 33 34 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. 36 If new_col_index[i]>-1 new_col_index[i] gives the column of i in 37 the returned submatrix. 38 */ 39 40 41 Paso_SparseMatrix* Paso_SparseMatrix_getSubmatrix(Paso_SparseMatrix* A,int n_row_sub,int n_col_sub, index_t* row_list,index_t* new_col_index){ 42 Paso_Pattern* sub_pattern=NULL; 43 Paso_SparseMatrix* out=NULL; 44 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0); 45 int i,k,tmp,m,subpattern_row; 46 int type=A->type; 47 Esys_resetError(); 48 if (A->type & MATRIX_FORMAT_CSC) { 49 Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_getSubmatrix: gathering submatrices supports CSR matrix format only."); 50 } else { 51 sub_pattern=Paso_Pattern_getSubpattern(A->pattern,n_row_sub,n_col_sub,row_list,new_col_index); 52 if (Esys_noError()) { 53 /* create the return object */ 54 out=Paso_SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size,TRUE); 55 if (Esys_noError()) { 56 #pragma omp parallel for private(i,k,m,subpattern_row,tmp) schedule(static) 57 for (i=0;ipattern->ptr[subpattern_row]-index_offset;kpattern->ptr[subpattern_row+1]-index_offset;++k) { 60 tmp=new_col_index[A->pattern->index[k]-index_offset]; 61 if (tmp>-1) { 62 #pragma ivdep 63 for (m=out->pattern->ptr[i]-index_offset;mpattern->ptr[i+1]-index_offset;++m) { 64 if (out->pattern->index[m]==tmp+index_offset) { 65 Paso_copyShortDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size])); 66 break; 67 } 68 } 69 } 70 } 71 } 72 } 73 } 74 Paso_Pattern_free(sub_pattern); 75 } 76 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;ipattern->ptr[i];iptrpattern->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;ipattern->ptr[i];iptrpattern->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;ipattern->ptr[i];iptrpattern->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;ipattern->ptr[i];iptrpattern->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;ipattern->ptr[i];iptrpattern->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;ipattern->ptr[i];iptrpattern->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 }