/[escript]/trunk/paso/src/Pattern_getSubpattern.c
ViewVC logotype

Contents of /trunk/paso/src/Pattern_getSubpattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1736 - (show annotations)
Fri Aug 29 02:23:16 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/plain
File size: 3087 byte(s)
This fixes a problem which is typically arising when using reduced order
with MPI and a "small" number of elements per processor. In this case it
can happen that the couple matrix is not using all entries sent to the
processor. The old implementations assumed that the indices will cover
the entire input. This assumption has been removed.


1
2 /* $Id: Pattern_getSubpattern.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: PatternPattern */
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 "Pattern.h"
29 #include "PasoUtil.h"
30
31 /**************************************************************/
32
33 /* creates Pattern */
34
35 Paso_Pattern* Paso_Pattern_getSubpattern(Paso_Pattern* pattern, \
36 int newNumRows, int newNumCols, index_t* row_list,index_t* new_col_index) {
37 index_t index_offset=(pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
38 Paso_Pattern*out=NULL;
39 index_t *ptr=NULL,*index=NULL,k,j,subpattern_row,tmp;
40 dim_t i;
41 Paso_resetError();
42
43 ptr=MEMALLOC(newNumRows+1,index_t);
44 if (! Paso_checkPtr(ptr)) {
45 #pragma omp parallel
46 {
47 #pragma omp for private(i) schedule(static)
48 for (i=0;i<newNumRows+1;++i) ptr[i]=0;
49
50 /* find the number column entries in each row */
51 #pragma omp for private(i,k,j,subpattern_row) schedule(static)
52 for (i=0;i<newNumRows;++i) {
53 j=0;
54 subpattern_row=row_list[i];
55 #pragma ivdep
56 for (k=pattern->ptr[subpattern_row]-index_offset;k<pattern->ptr[subpattern_row+1]-index_offset;++k)
57 if (new_col_index[pattern->index[k]-index_offset]>-1) j++;
58 ptr[i]=j;
59 }
60 }
61 /* accummulate ptr */
62 ptr[newNumRows]=Paso_Util_cumsum(newNumRows,ptr);
63 index=MEMALLOC(ptr[newNumRows],index_t);
64 if (Paso_checkPtr(index)) {
65 MEMFREE(ptr);
66 } else {
67 /* find the number column entries in each row */
68 #pragma omp parallel for private(i,k,j,subpattern_row,tmp) schedule(static)
69 for (i=0;i<newNumRows;++i) {
70 j=ptr[i];
71 subpattern_row=row_list[i];
72 #pragma ivdep
73 for (k=pattern->ptr[subpattern_row]-index_offset;k<pattern->ptr[subpattern_row+1]-index_offset;++k) {
74 tmp=new_col_index[pattern->index[k]-index_offset];
75 if (tmp>-1) {
76 index[j]=tmp;
77 ++j;
78 }
79 }
80 }
81 /* create return value */
82 out=Paso_Pattern_alloc(pattern->type,pattern->input_block_size,pattern->output_block_size,newNumRows,newNumCols,ptr,index);
83 if (! Paso_noError()) {
84 MEMFREE(index);
85 MEMFREE(ptr);
86 }
87 }
88 }
89 return out;
90 }

  ViewVC Help
Powered by ViewVC 1.1.26