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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 633 - (hide annotations)
Thu Mar 23 05:37:00 2006 UTC (13 years, 11 months ago) by dhawcroft
File MIME type: text/plain
File size: 3991 byte(s)


1 jgs 150 /* $Id$ */
2    
3 dhawcroft 631 /*
4     ********************************************************************************
5 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
6 dhawcroft 631 * *
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 jgs 150 /**************************************************************/
15    
16     /* Paso: SystemMatrixPatternPattern */
17    
18     /**************************************************************/
19    
20     /* Copyrights by ACcESS Australia 2003, 2004, 2005 */
21     /* Author: gross@access.edu.au */
22    
23     /**************************************************************/
24    
25     #include "Paso.h"
26     #include "SystemMatrixPattern.h"
27    
28     /**************************************************************/
29    
30     /* creates SystemMatrixPattern */
31    
32     Paso_SystemMatrixPattern* Paso_SystemMatrixPattern_unrollBlocks(Paso_SystemMatrixPattern* pattern, \
33 gross 415 int type, dim_t row_block_size,dim_t col_block_size) {
34     index_t index_offset_in=(pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
35     index_t index_offset_out=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
36    
37     if ((pattern->type & PATTERN_FORMAT_SYM) != (type & PATTERN_FORMAT_SYM)) {
38     Paso_setError(TYPE_ERROR,"Paso_SystemMatrixPattern_unrollBlocks: conversion between symmetric and non-symmetric is not implemented yet");
39     return NULL;
40     }
41 jgs 150 Paso_SystemMatrixPattern*out=NULL;
42     index_t *ptr=NULL,*index=NULL,iPtr;
43     dim_t i,j,k,l;
44     Paso_resetError();
45     dim_t block_size=row_block_size*col_block_size;
46     dim_t new_n_ptr=(pattern->n_ptr)*row_block_size;
47     dim_t new_len=(pattern->len)*block_size;
48    
49     ptr=MEMALLOC(new_n_ptr+1,index_t);
50     index=MEMALLOC(new_len,index_t);
51    
52    
53     if (! ( Paso_checkPtr(ptr) || Paso_checkPtr(index) ) ) {
54     #pragma omp parallel
55     {
56     #pragma omp for private(i) schedule(static)
57 gross 415 for (i=0;i<new_n_ptr+1;++i) ptr[i]=index_offset_out;
58 jgs 150
59     #pragma omp master
60 gross 415 ptr[new_n_ptr]=new_len+index_offset_out;
61 jgs 150
62     #pragma omp for private(i,k) schedule(static)
63     for (i=0;i<pattern->n_ptr;++i)
64 gross 415 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 jgs 150
66     #pragma omp for private(i,iPtr) schedule(static)
67     for (i=0;i<new_n_ptr;++i)
68 gross 432 for (iPtr=ptr[i]-index_offset_out;iPtr<ptr[i+1]-index_offset_out;++iPtr) index[iPtr]=index_offset_out;
69 jgs 150
70     #pragma omp for private(i,j,iPtr,k) schedule(static)
71     for (i=0;i<pattern->n_ptr;++i) {
72 gross 415 for (iPtr=pattern->ptr[i]-index_offset_in;iPtr<pattern->ptr[i+1]-index_offset_in;++iPtr) {
73 jgs 150 for (k=0;k<row_block_size;++k) {
74     for (j=0;j<col_block_size;++j) {
75 gross 425 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 jgs 150 }
77     }
78     }
79     }
80     }
81     /* create return value */
82 gross 415 out=Paso_SystemMatrixPattern_alloc(type,new_n_ptr,ptr,index);
83 jgs 150 }
84     if (! Paso_noError()) {
85     MEMFREE(index);
86     MEMFREE(ptr);
87     }
88     return out;
89     }
90     /*
91     * $Log$
92     * Revision 1.2 2005/09/15 03:44:39 jgs
93     * Merge of development branch dev-02 back to main trunk on 2005-09-15
94     *
95     * Revision 1.1.2.1 2005/09/05 06:29:47 gross
96     * These files have been extracted from finley to define a stand alone libray for iterative
97     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
98     * has not been tested yet.
99     *
100     *
101     */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26