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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3312 - (show annotations)
Tue Oct 26 07:54:58 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 3286 byte(s)
last step for a clean up version of the AMG
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 /* Author: Lutz Gross, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "SparseMatrix.h"
27
28
29 Paso_SparseMatrix* Paso_SparseMatrix_unroll(const Paso_SparseMatrixType type, const Paso_SparseMatrix* A) {
30
31 const dim_t row_block_size = A->row_block_size;
32 const dim_t col_block_size = A->col_block_size;
33 const dim_t block_size=A->block_size;
34 const dim_t n = A->numRows;
35 const index_t out_type = (type & MATRIX_FORMAT_BLK1) ? type : type + MATRIX_FORMAT_BLK1;
36 const index_t A_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
37 const index_t out_offset= (out_type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38 index_t iptr, *start_p, *where_p;
39 dim_t i, j, irb, icb, icol, irow, l_col, l_row;
40
41 Paso_SparseMatrix*out=NULL;
42
43 out = Paso_SparseMatrix_alloc(out_type, A->pattern, row_block_size, col_block_size, FALSE);
44
45 if (Esys_noError()) {
46 if (out->type & MATRIX_FORMAT_CSC) {
47
48 #pragma omp parallel for private(i, iptr, j, irb, irow, start_p, l_col, icb, icol, where_p) schedule(static)
49 for (i=0;i<n;++i) {
50 for (iptr=A->pattern->ptr[i]-A_offset; iptr<A->pattern->ptr[i+1]-A_offset; ++iptr) {
51 j=A->pattern->index[iptr]-A_offset;
52 for (icb=0; icb<col_block_size; ++icb) {
53 icol=j*col_block_size+icb;
54 start_p=&(out->pattern->index[out->pattern->ptr[icol]-out_offset]);
55 l_col=out->pattern->ptr[icol + 1]-out->pattern->ptr[icol];
56 for (irb=0; irb<row_block_size; ++irb) {
57 irow=row_block_size*i+irb+out_offset;
58 where_p=(index_t*)bsearch(&(irow), start_p, l_col,sizeof(index_t), Paso_comparIndex);
59 if (! (where_p == NULL) )
60 out->val[out->pattern->ptr[icol]-out_offset+(index_t)(where_p-start_p)] =
61 A->val[block_size*iptr+irb+row_block_size*icb];
62 }
63 }
64 }
65 }
66 } else {
67 #pragma omp parallel for private(i, iptr, j, irb, irow, start_p, l_row, icb, icol, where_p) schedule(static)
68 for (i=0;i<n;++i) {
69 for (iptr=A->pattern->ptr[i]-A_offset; iptr<A->pattern->ptr[i+1]-A_offset; ++iptr) {
70 j=A->pattern->index[iptr]-A_offset;
71 for (irb=0; irb<row_block_size; ++irb) {
72 irow=row_block_size*i+irb;
73 start_p=&(out->pattern->index[out->pattern->ptr[irow]-out_offset]);
74 l_row=out->pattern->ptr[irow + 1]-out->pattern->ptr[irow];
75 for (icb=0; icb<col_block_size; ++icb) {
76 icol=j*col_block_size+icb+out_offset;
77 where_p=(index_t*)bsearch(&(icol), start_p, l_row,sizeof(index_t), Paso_comparIndex);
78 if (! (where_p == NULL) )
79 out->val[out->pattern->ptr[irow]-out_offset+(index_t)(where_p-start_p)] =
80 A->val[block_size*iptr+irb+row_block_size*icb];
81 }
82 }
83 }
84 }
85 }
86 }
87 return out;
88 }

  ViewVC Help
Powered by ViewVC 1.1.26