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

Annotation of /trunk/paso/src/SparseMatrix_getTranspose.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4803 - (hide annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 10 months ago) by caltinay
File size: 4245 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

1 gross 3283
2 jfenwick 3981 /*****************************************************************************
3 gross 3283 *
4 jfenwick 4657 * Copyright (c) 2003-2014 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 gross 3283 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 gross 3283
17    
18 caltinay 4797 /****************************************************************************
19 gross 3283
20 gross 3303 Paso: transposed matrix
21 gross 3283
22 caltinay 4797 *****************************************************************************
23 gross 3283
24 caltinay 3642 Author: Lutz Gross, l.gross@uq.edu.au
25 gross 3283
26 caltinay 4797 *****************************************************************************/
27 gross 3283
28     #include "SparseMatrix.h"
29    
30 caltinay 4797 namespace paso {
31    
32     SparseMatrix* SparseMatrix_getTranspose(const SparseMatrix* A)
33 gross 3283 {
34    
35 caltinay 4803 Pattern *ATpattern=NULL;
36 caltinay 4797 SparseMatrix *AT=NULL;
37 gross 3283
38 gross 3303 const dim_t m=A->numCols;
39     const dim_t n=A->numRows;
40     const dim_t block_size=A->block_size;
41     const dim_t col_block_size_A = A->col_block_size;
42     const dim_t row_block_size_A = A->row_block_size;
43     register index_t iptr_AT, jptr_A, *start_p, *where_p, iptr2;
44     dim_t i;
45     register dim_t j, ib, irb, icb;
46 gross 3283
47 gross 3303 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(m);
48 gross 3283 for (i=0;i<n;++i) {
49 gross 3303 for (iptr2=A->pattern->ptr[i];iptr2<A->pattern->ptr[i+1]; ++iptr2) {
50     j=A->pattern->index[iptr2];
51 gross 3283 Paso_IndexListArray_insertIndex(index_list,j,i);
52 gross 3303 }
53     }
54    
55 caltinay 4803 ATpattern=Pattern_fromIndexListArray(0,index_list,0,n,0);
56 gross 3303 Paso_IndexListArray_free(index_list);
57 caltinay 4797 AT=SparseMatrix_alloc(A->type,ATpattern,col_block_size_A,row_block_size_A,FALSE);
58 caltinay 4803 Pattern_free(ATpattern);
59 gross 3303
60     if ( ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) && (block_size == 1 ) ) ||
61     ( (row_block_size_A == 1 ) && (col_block_size_A == 1) ) ) {
62     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p)
63     for (i=0;i<m;++i) {
64     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
65     j=AT->pattern->index[iptr_AT];
66     jptr_A=A->pattern->ptr[j];
67     start_p=&(A->pattern->index[jptr_A]);
68     where_p=(index_t*)bsearch(&i, start_p,
69     A->pattern->ptr[j + 1]-jptr_A,
70     sizeof(index_t),
71 caltinay 4803 comparIndex);
72 gross 3303 if (! (where_p == NULL) ) { /* this should always be the case */
73     jptr_A += (index_t)(where_p-start_p);
74     AT->val[iptr_AT]=A->val[jptr_A];
75     }
76     }
77 gross 3283 }
78 gross 3303 } else {
79     if (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
80     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, ib)
81     for (i=0;i<m;++i) {
82     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
83     j=AT->pattern->index[iptr_AT];
84     jptr_A=A->pattern->ptr[j];
85     start_p=&(A->pattern->index[jptr_A]);
86     where_p=(index_t*)bsearch(&i, start_p,
87     A->pattern->ptr[j + 1]-jptr_A,
88     sizeof(index_t),
89 caltinay 4803 comparIndex);
90 gross 3303 if (! (where_p == NULL) ) { /* this should always be the case */
91     jptr_A += (index_t)(where_p-start_p);
92     for (ib=0; ib < block_size; ++ib ) AT->val[iptr_AT*block_size+ib]=A->val[jptr_A*block_size+ib];
93     }
94     }
95 gross 3283 }
96 gross 3303 } else {
97     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, irb, icb)
98     for (i=0;i<m;++i) {
99     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
100     j=AT->pattern->index[iptr_AT];
101     jptr_A=A->pattern->ptr[j];
102     start_p=&(A->pattern->index[jptr_A]);
103     where_p=(index_t*)bsearch(&i, start_p,
104     A->pattern->ptr[j + 1]-jptr_A,
105     sizeof(index_t),
106 caltinay 4803 comparIndex);
107 gross 3303 if (! (where_p == NULL) ) { /* this should always be the case */
108     jptr_A += (index_t)(where_p-start_p);
109     for (irb =0 ; irb < row_block_size_A; ++irb) {
110     for (icb =0 ; icb < col_block_size_A; ++icb) {
111     AT->val[iptr_AT*block_size+icb+col_block_size_A*irb]=A->val[jptr_A*block_size+irb+row_block_size_A*icb];
112     }
113     }
114     }
115     }
116     }
117     }
118     }
119     return AT;
120 caltinay 3642 }
121 caltinay 4797
122     } // namespace paso
123    

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/SparseMatrix_getTranspose.cpp:3531-3826 /branches/lapack2681/paso/src/SparseMatrix_getTranspose.cpp:2682-2741 /branches/pasowrap/paso/src/SparseMatrix_getTranspose.cpp:3661-3674 /branches/py3_attempt2/paso/src/SparseMatrix_getTranspose.cpp:3871-3891 /branches/restext/paso/src/SparseMatrix_getTranspose.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/SparseMatrix_getTranspose.cpp:3669-3791 /branches/stage3.0/paso/src/SparseMatrix_getTranspose.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/SparseMatrix_getTranspose.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/SparseMatrix_getTranspose.cpp:3517-3974 /release/3.0/paso/src/SparseMatrix_getTranspose.cpp:2591-2601 /trunk/paso/src/SparseMatrix_getTranspose.cpp:4257-4344 /trunk/ripley/test/python/paso/src/SparseMatrix_getTranspose.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26