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

Diff of /trunk/paso/src/SystemMatrixPattern_unrollBlocks.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    /*
4    ********************************************************************************
5    *               Copyright   2006 by ACcESS MNRF                                *
6    *                                                                              *
7    *                 http://www.access.edu.au                                     *
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  /* Paso: SystemMatrixPatternPattern */  /* Paso: SystemMatrixPatternPattern */
# Line 19  Line 30 
30  /* creates SystemMatrixPattern  */  /* creates SystemMatrixPattern  */
31    
32  Paso_SystemMatrixPattern* Paso_SystemMatrixPattern_unrollBlocks(Paso_SystemMatrixPattern* pattern, \  Paso_SystemMatrixPattern* Paso_SystemMatrixPattern_unrollBlocks(Paso_SystemMatrixPattern* pattern, \
33                                             dim_t row_block_size,dim_t col_block_size) {                                             int type, dim_t row_block_size,dim_t col_block_size) {
34    Paso_SystemMatrixPattern*out=NULL;    Paso_SystemMatrixPattern*out=NULL;
35    index_t *ptr=NULL,*index=NULL,iPtr;    index_t *ptr=NULL,*index=NULL,iPtr;
36    dim_t i,j,k,l;    dim_t i,j,k,l, block_size, new_n_ptr, new_len;
37      index_t index_offset_in=(pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
38      index_t index_offset_out=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
39      
40      if ((pattern->type & PATTERN_FORMAT_SYM) != (type & PATTERN_FORMAT_SYM)) {
41          Paso_setError(TYPE_ERROR,"Paso_SystemMatrixPattern_unrollBlocks: conversion between symmetric and non-symmetric is not implemented yet");
42          return NULL;
43      }
44    Paso_resetError();    Paso_resetError();
45    dim_t block_size=row_block_size*col_block_size;    block_size=row_block_size*col_block_size;
46    dim_t new_n_ptr=(pattern->n_ptr)*row_block_size;    new_n_ptr=(pattern->n_ptr)*row_block_size;
47    dim_t new_len=(pattern->len)*block_size;    new_len=(pattern->len)*block_size;
48    
49    ptr=MEMALLOC(new_n_ptr+1,index_t);    ptr=MEMALLOC(new_n_ptr+1,index_t);
50    index=MEMALLOC(new_len,index_t);    index=MEMALLOC(new_len,index_t);
# Line 36  Paso_SystemMatrixPattern* Paso_SystemMat Line 54  Paso_SystemMatrixPattern* Paso_SystemMat
54       #pragma omp parallel       #pragma omp parallel
55       {       {
56          #pragma omp for private(i) schedule(static)          #pragma omp for private(i) schedule(static)
57          for (i=0;i<new_n_ptr+1;++i) ptr[i]=0;          for (i=0;i<new_n_ptr+1;++i) ptr[i]=index_offset_out;
58    
59          #pragma omp master          #pragma omp master
60          ptr[new_n_ptr]=new_len;          ptr[new_n_ptr]=new_len+index_offset_out;
61    
62          #pragma omp for private(i,k) schedule(static)          #pragma omp for private(i,k) schedule(static)
63          for (i=0;i<pattern->n_ptr;++i)          for (i=0;i<pattern->n_ptr;++i)
64              for (k=0;k<row_block_size;++k) ptr[i*row_block_size+k]=(pattern->ptr[i]-PTR_OFFSET)*block_size+(pattern->ptr[i+1]-pattern->ptr[i])*col_block_size*k;              for (k=0;k<row_block_size;++k) ptr[i*row_block_size+k]=(pattern->ptr[i]-index_offset_in)*block_size+(pattern->ptr[i+1]-pattern->ptr[i])*col_block_size*k+index_offset_out;
65                        
66          #pragma omp for private(i,iPtr) schedule(static)          #pragma omp for private(i,iPtr) schedule(static)
67          for (i=0;i<new_n_ptr;++i)          for (i=0;i<new_n_ptr;++i)
68              for (iPtr=ptr[i];iPtr<ptr[i+1];++iPtr) index[iPtr]=0;              for (iPtr=ptr[i]-index_offset_out;iPtr<ptr[i+1]-index_offset_out;++iPtr) index[iPtr]=index_offset_out;
69    
70          #pragma omp for private(i,j,iPtr,k) schedule(static)          #pragma omp for private(i,j,iPtr,k) schedule(static)
71          for (i=0;i<pattern->n_ptr;++i) {          for (i=0;i<pattern->n_ptr;++i) {
72             for (iPtr=pattern->ptr[i];iPtr<pattern->ptr[i+1];++iPtr)  {             for (iPtr=pattern->ptr[i]-index_offset_in;iPtr<pattern->ptr[i+1]-index_offset_in;++iPtr)  {
73                for (k=0;k<row_block_size;++k) {                for (k=0;k<row_block_size;++k) {
74                   for (j=0;j<col_block_size;++j) {                   for (j=0;j<col_block_size;++j) {
75                      index[ptr[i*row_block_size+k]+(iPtr-pattern->ptr[i])*col_block_size+j]=(pattern->index[iPtr]-INDEX_OFFSET)*col_block_size+j;                      index[ptr[i*row_block_size+k]-index_offset_out+(iPtr-(pattern->ptr[i]-index_offset_in))*col_block_size+j]=(pattern->index[iPtr]-index_offset_in)*col_block_size+j+index_offset_out;
76                   }                   }
77                }                }
78             }             }
79          }          }
80       }       }
81       /* create return value */       /* create return value */
82       out=Paso_SystemMatrixPattern_alloc(new_n_ptr,ptr,index);       out=Paso_SystemMatrixPattern_alloc(type,new_n_ptr,ptr,index);
83    }      }  
84    if (! Paso_noError()) {    if (! Paso_noError()) {
85       MEMFREE(index);       MEMFREE(index);

Legend:
Removed from v.155  
changed lines
  Added in v.1028

  ViewVC Help
Powered by ViewVC 1.1.26