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

Revision 1315 - (hide annotations)
Tue Sep 25 02:41:13 2007 UTC (13 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 3076 byte(s)
Copied more files from MPI branch to trunk

 1 ksteube 1315 2 /* \$Id: SparseMatrix_getSubmatrix.c 1306 2007-09-18 05:51:09Z ksteube \$ */ 3 4 /******************************************************* 5 * 6 * Copyright 2003-2007 by ACceSS MNRF 7 * Copyright 2007 by University of Queensland 8 * 9 * http://esscc.uq.edu.au 10 * Primary Business: Queensland, Australia 11 * Licensed under the Open Software License version 3.0 12 * http://www.opensource.org/licenses/osl-3.0.php 13 * 14 *******************************************************/ 15 16 /**************************************************************/ 17 18 /* Paso: SparseMatrix */ 19 20 /**************************************************************/ 21 22 /* Copyrights by ACcESS Australia 2003, 2004,2005 */ 23 /* Author: gross@access.edu.au */ 24 25 /**************************************************************/ 26 27 #include "Paso.h" 28 #include "SparseMatrix.h" 29 #include "PasoUtil.h" 30 31 /************************************************************** 32 33 returns the submatrix of A where rows are gathered by index row_list 34 and columns are selected by non-negative values of new_col_index. 35 if new_col_index[i]>-1 new_col_index[i] gives the column of i in 36 the returned submatrix 37 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 Paso_resetError(); 48 if (A->type & MATRIX_FORMAT_CSC) { 49 Paso_setError(TYPE_ERROR,"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 (Paso_noError()) { 53 /* create the return object */ 54 out=Paso_SparseMatrix_alloc(type,sub_pattern,A->row_block_size,A->col_block_size); 55 if (Paso_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_copyDouble(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 }