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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26