# Contents of /branches/doubleplusgood/paso/src/SparseMatrix_getSubmatrix.cpp

Revision 2666 - (show annotations)
Thu Sep 17 00:10:00 2009 UTC (10 years, 2 months ago) by artak
Original Path: trunk/paso/src/SparseMatrix_getSubmatrix.c
File MIME type: text/plain
File size: 5985 byte(s)
```OpenMP bug fixed and more openMP code added for cutting blocks.
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2009 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 14 15 /**************************************************************/ 16 17 /* Paso: SparseMatrix */ 18 19 /**************************************************************/ 20 21 /* Copyrights by ACcESS Australia 2003, 2004,2005 */ 22 /* Author: Lutz Gross, l.gross@uq.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,TRUE); 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 Paso_copyShortDouble(A->block_size,&(A->val[k*A->block_size]),&(out->val[m*A->block_size])); 65 break; 66 } 67 } 68 } 69 } 70 } 71 } 72 } 73 Paso_Pattern_free(sub_pattern); 74 } 75 return out; 76 } 77 78 Paso_SparseMatrix* Paso_SparseMatrix_getBlock(Paso_SparseMatrix* A, int blockid){ 79 80 dim_t blocksize=A->row_block_size,i; 81 dim_t n=A->numRows; 82 index_t iptr; 83 84 Paso_SparseMatrix* out=NULL; 85 86 out=Paso_SparseMatrix_alloc(A->type,A->pattern, 1, 1, 0); 87 88 if (blocksize==1) { 89 if (blockid==1) { 90 #pragma omp parallel for private(i,iptr) schedule(static) 91 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 93 out->val[iptr]=A->val[iptr]; 94 } 95 } 96 } 97 else { 98 Paso_setError(VALUE_ERROR,"Requested and actual block sizes do not match."); 99 } 100 } 101 else if (blocksize==2) { 102 if (blockid==1) { 103 #pragma omp parallel for private(i,iptr) schedule(static) 104 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 106 out->val[iptr]=A->val[4*iptr]; 107 } 108 } 109 } 110 else if (blockid==2) { 111 #pragma omp parallel for private(i,iptr) schedule(static) 112 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 114 out->val[iptr]=A->val[4*iptr+3]; 115 } 116 } 117 } 118 else { 119 Paso_setError(VALUE_ERROR,"Requested and actual block sizes do not match."); 120 } 121 } 122 else if (blocksize==3) { 123 if (blockid==1) { 124 #pragma omp parallel for private(i,iptr) schedule(static) 125 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 127 out->val[iptr]=A->val[9*iptr]; 128 } 129 } 130 } 131 else if (blockid==2) { 132 #pragma omp parallel for private(i,iptr) schedule(static) 133 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 135 out->val[iptr]=A->val[9*iptr+4]; 136 } 137 } 138 } 139 else if (blockid==3) { 140 #pragma omp parallel for private(i,iptr) schedule(static) 141 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 143 out->val[iptr]=A->val[9*iptr+8]; 144 } 145 } 146 } 147 else { 148 Paso_setError(VALUE_ERROR,"Requested and actual block sizes do not match."); 149 } 150 } 151 152 return out; 153 }