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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (9 years, 4 months ago) by gross
Original Path: trunk/paso/src/SparseMatrix_getTranspose.c
File MIME type: text/plain
File size: 4053 byte(s)
AMG reengineered, BUG is Smoother fixed.


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: inverts the main diagonal of a SparseMatrix: */
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* P)
29 {
30
31 Paso_Pattern *outpattern=NULL;
32 Paso_SparseMatrix *out=NULL;
33
34 const dim_t C=P->numCols;
35 const dim_t F=P->numRows-C;
36 const dim_t n=C+F;
37 const dim_t block_size=P->row_block_size;
38 dim_t i,j,k=0;
39 index_t iptr,jptr;
40 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(C);
41
42
43 for (i=0;i<n;++i) {
44 for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {
45 j=P->pattern->index[iptr];
46 Paso_IndexListArray_insertIndex(index_list,j,i);
47 }
48 }
49
50 outpattern=Paso_Pattern_fromIndexListArray(0,index_list,0,C+F,0);
51 out=Paso_SparseMatrix_alloc(P->type,outpattern,block_size,block_size,FALSE);
52
53
54 if (block_size==1) {
55 for (i=0;i<out->numRows;++i) {
56 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
57 j=out->pattern->index[iptr];
58 /*This can be replaced by bsearch!!*/
59 for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
60 k=P->pattern->index[jptr];
61 if(k==i) {
62 out->val[iptr]=P->val[jptr];
63 }
64 }
65 }
66 }
67 } else if (block_size==2) {
68 for (i=0;i<out->numRows;++i) {
69 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
70 j=out->pattern->index[iptr];
71 /*This can be replaced by bsearch!!*/
72 for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
73 k=P->pattern->index[jptr];
74 if(k==i) {
75 out->val[iptr*block_size*block_size]=P->val[jptr*block_size*block_size];
76 out->val[iptr*block_size*block_size+1]=P->val[jptr*block_size*block_size+2];
77 out->val[iptr*block_size*block_size+2]=P->val[jptr*block_size*block_size+1];
78 out->val[iptr*block_size*block_size+3]=P->val[jptr*block_size*block_size+3];
79 }
80 }
81 }
82 }
83 } else if (block_size==3) {
84 for (i=0;i<out->numRows;++i) {
85 for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
86 j=out->pattern->index[iptr];
87 /*This can be replaced by bsearch!!*/
88 for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
89 k=P->pattern->index[jptr];
90 if(k==i) {
91 out->val[iptr*block_size*block_size]=P->val[jptr*block_size*block_size];
92 out->val[iptr*block_size*block_size+1]=P->val[jptr*block_size*block_size+3];
93 out->val[iptr*block_size*block_size+2]=P->val[jptr*block_size*block_size+6];
94 out->val[iptr*block_size*block_size+3]=P->val[jptr*block_size*block_size+1];
95 out->val[iptr*block_size*block_size+4]=P->val[jptr*block_size*block_size+4];
96 out->val[iptr*block_size*block_size+5]=P->val[jptr*block_size*block_size+7];
97 out->val[iptr*block_size*block_size+6]=P->val[jptr*block_size*block_size+2];
98 out->val[iptr*block_size*block_size+7]=P->val[jptr*block_size*block_size+5];
99 out->val[iptr*block_size*block_size+8]=P->val[jptr*block_size*block_size+8];
100 }
101 }
102 }
103 }
104 }
105
106 /* clean up */
107 Paso_IndexListArray_free(index_list);
108 Paso_Pattern_free(outpattern);
109 return out;
110 }
111

  ViewVC Help
Powered by ViewVC 1.1.26