/[escript]/branches/doubleplusgood/paso/src/SparseMatrix_unroll.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SparseMatrix_unroll.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4261 - (show annotations)
Wed Feb 27 06:09:33 2013 UTC (6 years, 9 months ago) by jfenwick
File size: 3472 byte(s)
Initial all c++ build.
But ... there are now reinterpret_cast<>'s
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: SparseMatrix */
20
21 /************************************************************************************/
22
23 /* Author: Lutz Gross, l.gross@uq.edu.au */
24
25 /************************************************************************************/
26
27 #include "Paso.h"
28 #include "SparseMatrix.h"
29
30
31 Paso_SparseMatrix* Paso_SparseMatrix_unroll(const Paso_SparseMatrixType type, const Paso_SparseMatrix* A) {
32
33 const dim_t row_block_size = A->row_block_size;
34 const dim_t col_block_size = A->col_block_size;
35 const dim_t block_size=A->block_size;
36 const dim_t n = A->numRows;
37 const index_t out_type = (type & MATRIX_FORMAT_BLK1) ? type : type + MATRIX_FORMAT_BLK1;
38 const index_t A_offset=(A->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
39 const index_t out_offset= (out_type & MATRIX_FORMAT_OFFSET1 ? 1:0);
40 index_t iptr, *start_p, *where_p;
41 dim_t i, j, irb, icb, icol, irow, l_col, l_row;
42
43 Paso_SparseMatrix*out=NULL;
44
45 out = Paso_SparseMatrix_alloc(out_type, A->pattern, row_block_size, col_block_size, FALSE);
46
47 if (Esys_noError()) {
48 if (out->type & MATRIX_FORMAT_CSC) {
49
50 #pragma omp parallel for private(i, iptr, j, irb, irow, start_p, l_col, icb, icol, where_p) schedule(static)
51 for (i=0;i<n;++i) {
52 for (iptr=A->pattern->ptr[i]-A_offset; iptr<A->pattern->ptr[i+1]-A_offset; ++iptr) {
53 j=A->pattern->index[iptr]-A_offset;
54 for (icb=0; icb<col_block_size; ++icb) {
55 icol=j*col_block_size+icb;
56 start_p=&(out->pattern->index[out->pattern->ptr[icol]-out_offset]);
57 l_col=out->pattern->ptr[icol + 1]-out->pattern->ptr[icol];
58 for (irb=0; irb<row_block_size; ++irb) {
59 irow=row_block_size*i+irb+out_offset;
60 where_p=(index_t*)bsearch(&(irow), start_p, l_col,sizeof(index_t), Paso_comparIndex);
61 if (! (where_p == NULL) )
62 out->val[out->pattern->ptr[icol]-out_offset+(index_t)(where_p-start_p)] =
63 A->val[block_size*iptr+irb+row_block_size*icb];
64 }
65 }
66 }
67 }
68 } else {
69 #pragma omp parallel for private(i, iptr, j, irb, irow, start_p, l_row, icb, icol, where_p) schedule(static)
70 for (i=0;i<n;++i) {
71 for (iptr=A->pattern->ptr[i]-A_offset; iptr<A->pattern->ptr[i+1]-A_offset; ++iptr) {
72 j=A->pattern->index[iptr]-A_offset;
73 for (irb=0; irb<row_block_size; ++irb) {
74 irow=row_block_size*i+irb;
75 start_p=&(out->pattern->index[out->pattern->ptr[irow]-out_offset]);
76 l_row=out->pattern->ptr[irow + 1]-out->pattern->ptr[irow];
77 for (icb=0; icb<col_block_size; ++icb) {
78 icol=j*col_block_size+icb+out_offset;
79 where_p=(index_t*)bsearch(&(icol), start_p, l_row,sizeof(index_t), Paso_comparIndex);
80 if (! (where_p == NULL) )
81 out->val[out->pattern->ptr[irow]-out_offset+(index_t)(where_p-start_p)] =
82 A->val[block_size*iptr+irb+row_block_size*icb];
83 }
84 }
85 }
86 }
87 }
88 }
89 return out;
90 }

  ViewVC Help
Powered by ViewVC 1.1.26