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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4261 - (show annotations)
Wed Feb 27 06:09:33 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 4208 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: transposed matrix
20
21 ************************************************************************************
22
23 Author: Lutz Gross, l.gross@uq.edu.au
24
25 ************************************************************************************/
26
27 #include "SparseMatrix.h"
28
29 Paso_SparseMatrix* Paso_SparseMatrix_getTranspose(Paso_SparseMatrix* A)
30 {
31
32 Paso_Pattern *ATpattern=NULL;
33 Paso_SparseMatrix *AT=NULL;
34
35 const dim_t m=A->numCols;
36 const dim_t n=A->numRows;
37 const dim_t block_size=A->block_size;
38 const dim_t col_block_size_A = A->col_block_size;
39 const dim_t row_block_size_A = A->row_block_size;
40 register index_t iptr_AT, jptr_A, *start_p, *where_p, iptr2;
41 dim_t i;
42 register dim_t j, ib, irb, icb;
43
44 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(m);
45 for (i=0;i<n;++i) {
46 for (iptr2=A->pattern->ptr[i];iptr2<A->pattern->ptr[i+1]; ++iptr2) {
47 j=A->pattern->index[iptr2];
48 Paso_IndexListArray_insertIndex(index_list,j,i);
49 }
50 }
51
52 ATpattern=Paso_Pattern_fromIndexListArray(0,index_list,0,n,0);
53 Paso_IndexListArray_free(index_list);
54 AT=Paso_SparseMatrix_alloc(A->type,ATpattern,col_block_size_A,row_block_size_A,FALSE);
55 Paso_Pattern_free(ATpattern);
56
57 if ( ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) && (block_size == 1 ) ) ||
58 ( (row_block_size_A == 1 ) && (col_block_size_A == 1) ) ) {
59 #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p)
60 for (i=0;i<m;++i) {
61 for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
62 j=AT->pattern->index[iptr_AT];
63 jptr_A=A->pattern->ptr[j];
64 start_p=&(A->pattern->index[jptr_A]);
65 where_p=(index_t*)bsearch(&i, start_p,
66 A->pattern->ptr[j + 1]-jptr_A,
67 sizeof(index_t),
68 Paso_comparIndex);
69 if (! (where_p == NULL) ) { /* this should always be the case */
70 jptr_A += (index_t)(where_p-start_p);
71 AT->val[iptr_AT]=A->val[jptr_A];
72 }
73 }
74 }
75 } else {
76 if (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
77 #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, ib)
78 for (i=0;i<m;++i) {
79 for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
80 j=AT->pattern->index[iptr_AT];
81 jptr_A=A->pattern->ptr[j];
82 start_p=&(A->pattern->index[jptr_A]);
83 where_p=(index_t*)bsearch(&i, start_p,
84 A->pattern->ptr[j + 1]-jptr_A,
85 sizeof(index_t),
86 Paso_comparIndex);
87 if (! (where_p == NULL) ) { /* this should always be the case */
88 jptr_A += (index_t)(where_p-start_p);
89 for (ib=0; ib < block_size; ++ib ) AT->val[iptr_AT*block_size+ib]=A->val[jptr_A*block_size+ib];
90 }
91 }
92 }
93 } else {
94 #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, irb, icb)
95 for (i=0;i<m;++i) {
96 for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
97 j=AT->pattern->index[iptr_AT];
98 jptr_A=A->pattern->ptr[j];
99 start_p=&(A->pattern->index[jptr_A]);
100 where_p=(index_t*)bsearch(&i, start_p,
101 A->pattern->ptr[j + 1]-jptr_A,
102 sizeof(index_t),
103 Paso_comparIndex);
104 if (! (where_p == NULL) ) { /* this should always be the case */
105 jptr_A += (index_t)(where_p-start_p);
106 for (irb =0 ; irb < row_block_size_A; ++irb) {
107 for (icb =0 ; icb < col_block_size_A; ++icb) {
108 AT->val[iptr_AT*block_size+icb+col_block_size_A*irb]=A->val[jptr_A*block_size+irb+row_block_size_A*icb];
109 }
110 }
111 }
112 }
113 }
114 }
115 }
116 return AT;
117 }

  ViewVC Help
Powered by ViewVC 1.1.26