/[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 4154 - (hide annotations)
Tue Jan 22 09:30:23 2013 UTC (7 years, 1 month ago) by jfenwick
Original Path: trunk/paso/src/SparseMatrix_getTranspose.c
File MIME type: text/plain
File size: 4248 byte(s)
Round 1 of copyright fixes
1 gross 3283
2 jfenwick 3981 /*****************************************************************************
3 gross 3283 *
4 jfenwick 4154 * Copyright (c) 2003-2013 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     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 gross 3283
16    
17 jfenwick 3981 /************************************************************************************
18 gross 3283
19 gross 3303 Paso: transposed matrix
20 gross 3283
21 jfenwick 3981 ************************************************************************************
22 gross 3283
23 gross 3303 Copyrights by ACcESS Australia 2010
24 caltinay 3642 Author: Lutz Gross, l.gross@uq.edu.au
25 gross 3283
26 jfenwick 3981 ************************************************************************************/
27 gross 3283
28     #include "SparseMatrix.h"
29    
30 gross 3303 Paso_SparseMatrix* Paso_SparseMatrix_getTranspose(Paso_SparseMatrix* A)
31 gross 3283 {
32    
33 gross 3303 Paso_Pattern *ATpattern=NULL;
34     Paso_SparseMatrix *AT=NULL;
35 gross 3283
36 gross 3303 const dim_t m=A->numCols;
37     const dim_t n=A->numRows;
38     const dim_t block_size=A->block_size;
39     const dim_t col_block_size_A = A->col_block_size;
40     const dim_t row_block_size_A = A->row_block_size;
41     register index_t iptr_AT, jptr_A, *start_p, *where_p, iptr2;
42     dim_t i;
43     register dim_t j, ib, irb, icb;
44 gross 3283
45 gross 3303 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(m);
46 gross 3283 for (i=0;i<n;++i) {
47 gross 3303 for (iptr2=A->pattern->ptr[i];iptr2<A->pattern->ptr[i+1]; ++iptr2) {
48     j=A->pattern->index[iptr2];
49 gross 3283 Paso_IndexListArray_insertIndex(index_list,j,i);
50 gross 3303 }
51     }
52    
53     ATpattern=Paso_Pattern_fromIndexListArray(0,index_list,0,n,0);
54     Paso_IndexListArray_free(index_list);
55     AT=Paso_SparseMatrix_alloc(A->type,ATpattern,col_block_size_A,row_block_size_A,FALSE);
56     Paso_Pattern_free(ATpattern);
57    
58     if ( ( (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) && (block_size == 1 ) ) ||
59     ( (row_block_size_A == 1 ) && (col_block_size_A == 1) ) ) {
60     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p)
61     for (i=0;i<m;++i) {
62     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
63     j=AT->pattern->index[iptr_AT];
64     jptr_A=A->pattern->ptr[j];
65     start_p=&(A->pattern->index[jptr_A]);
66     where_p=(index_t*)bsearch(&i, start_p,
67     A->pattern->ptr[j + 1]-jptr_A,
68     sizeof(index_t),
69     Paso_comparIndex);
70     if (! (where_p == NULL) ) { /* this should always be the case */
71     jptr_A += (index_t)(where_p-start_p);
72     AT->val[iptr_AT]=A->val[jptr_A];
73     }
74     }
75 gross 3283 }
76 gross 3303 } else {
77     if (A->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
78     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, ib)
79     for (i=0;i<m;++i) {
80     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
81     j=AT->pattern->index[iptr_AT];
82     jptr_A=A->pattern->ptr[j];
83     start_p=&(A->pattern->index[jptr_A]);
84     where_p=(index_t*)bsearch(&i, start_p,
85     A->pattern->ptr[j + 1]-jptr_A,
86     sizeof(index_t),
87     Paso_comparIndex);
88     if (! (where_p == NULL) ) { /* this should always be the case */
89     jptr_A += (index_t)(where_p-start_p);
90     for (ib=0; ib < block_size; ++ib ) AT->val[iptr_AT*block_size+ib]=A->val[jptr_A*block_size+ib];
91     }
92     }
93 gross 3283 }
94 gross 3303 } else {
95     #pragma omp parallel for private(i, iptr_AT, j, jptr_A, start_p, where_p, irb, icb)
96     for (i=0;i<m;++i) {
97     for (iptr_AT=AT->pattern->ptr[i];iptr_AT<AT->pattern->ptr[i+1]; ++iptr_AT) {
98     j=AT->pattern->index[iptr_AT];
99     jptr_A=A->pattern->ptr[j];
100     start_p=&(A->pattern->index[jptr_A]);
101     where_p=(index_t*)bsearch(&i, start_p,
102     A->pattern->ptr[j + 1]-jptr_A,
103     sizeof(index_t),
104     Paso_comparIndex);
105     if (! (where_p == NULL) ) { /* this should always be the case */
106     jptr_A += (index_t)(where_p-start_p);
107     for (irb =0 ; irb < row_block_size_A; ++irb) {
108     for (icb =0 ; icb < col_block_size_A; ++icb) {
109     AT->val[iptr_AT*block_size+icb+col_block_size_A*irb]=A->val[jptr_A*block_size+irb+row_block_size_A*icb];
110     }
111     }
112     }
113     }
114     }
115     }
116     }
117     return AT;
118 caltinay 3642 }

  ViewVC Help
Powered by ViewVC 1.1.26