Annotation of /trunk/paso/src/SparseMatrix_getSubmatrix.c

Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (12 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 2982 byte(s)
```Copyright updated in all files

```
 1 ksteube 1315 2 /******************************************************* 3 ksteube 1811 * 4 * Copyright (c) 2003-2008 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 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 ksteube 1315 14 ksteube 1811 15 ksteube 1315 /**************************************************************/ 16 17 /* Paso: SparseMatrix */ 18 19 /**************************************************************/ 20 21 /* Copyrights by ACcESS Australia 2003, 2004,2005 */ 22 /* Author: gross@access.edu.au */ 23 24 /**************************************************************/ 25 26 #include "Paso.h" 27 #include "SparseMatrix.h" 28 #include "PasoUtil.h" 29 30 /************************************************************** 31 32 returns the submatrix of A where rows are gathered by index row_list 33 and columns are selected by non-negative values of new_col_index. 34 if new_col_index[i]>-1 new_col_index[i] gives the column of i in 35 the returned submatrix 36 37 */ 38 39 40 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){ 41 Paso_Pattern* sub_pattern=NULL; 42 Paso_SparseMatrix* out=NULL; 43 index_t index_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0); 44 int i,k,tmp,m,subpattern_row; 45 int type=A->type; 46 Paso_resetError(); 47 if (A->type & MATRIX_FORMAT_CSC) { 48 Paso_setError(TYPE_ERROR,"gathering submatrices supports CSR matrix format only."); 49 } else { 50 sub_pattern=Paso_Pattern_getSubpattern(A->pattern,n_row_sub,n_col_sub,row_list,new_col_index); 51 if (Paso_noError()) { 52 /* create the return object */ 53 out=Paso_SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size); 54 if (Paso_noError()) { 55 #pragma omp parallel for private(i,k,m,subpattern_row,tmp) schedule(static) 56 for (i=0;ipattern->ptr[subpattern_row]-index_offset;kpattern->ptr[subpattern_row+1]-index_offset;++k) { 59 tmp=new_col_index[A->pattern->index[k]-index_offset]; 60 if (tmp>-1) { 61 #pragma ivdep 62 for (m=out->pattern->ptr[i]-index_offset;mpattern->ptr[i+1]-index_offset;++m) { 63 if (out->pattern->index[m]==tmp+index_offset) { 64 gross 1639 Paso_copyShortDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size])); 65 ksteube 1315 break; 66 } 67 } 68 } 69 } 70 } 71 } 72 } 73 Paso_Pattern_free(sub_pattern); 74 } 75 return out; 76 }