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

Diff of /branches/doubleplusgood/paso/src/SparseMatrix_getTranspose.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3283 by gross, Mon Oct 18 22:39:28 2010 UTC revision 3303 by gross, Mon Oct 25 04:33:31 2010 UTC
# Line 12  Line 12 
12  *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************
16    
17  /* Paso: inverts the main diagonal of a SparseMatrix: */    Paso: transposed matrix
18    
19  /**************************************************************/  **************************************************************
20    
21  /* Copyrights by ACcESS Australia 2010 */     Copyrights by ACcESS Australia 2010
22  /* Author: Lutz Gross, l.gross@uq.edu.au */   Author: Lutz Gross, l.gross@uq.edu.au
23    
24  /**************************************************************/  **************************************************************/
25    
26  #include "SparseMatrix.h"  #include "SparseMatrix.h"
27    
28  Paso_SparseMatrix* Paso_SparseMatrix_getTranspose(Paso_SparseMatrix* P)  Paso_SparseMatrix* Paso_SparseMatrix_getTranspose(Paso_SparseMatrix* A)
29  {  {
30        
31     Paso_Pattern *outpattern=NULL;     Paso_Pattern *ATpattern=NULL;
32     Paso_SparseMatrix *out=NULL;     Paso_SparseMatrix *AT=NULL;
     
    const dim_t C=P->numCols;  
    const dim_t F=P->numRows-C;  
    const dim_t n=C+F;  
    const dim_t block_size=P->row_block_size;  
    dim_t i,j,k=0;  
    index_t iptr,jptr;  
    Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(C);  
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) {     for (i=0;i<n;++i) {
45        for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {        for (iptr2=A->pattern->ptr[i];iptr2<A->pattern->ptr[i+1]; ++iptr2) {
46       j=P->pattern->index[iptr];       j=A->pattern->index[iptr2];
47       Paso_IndexListArray_insertIndex(index_list,j,i);       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       }       }
       
      outpattern=Paso_Pattern_fromIndexListArray(0,index_list,0,C+F,0);  
      out=Paso_SparseMatrix_alloc(P->type,outpattern,block_size,block_size,FALSE);  
       
       
      if (block_size==1) {  
         for (i=0;i<out->numRows;++i) {  
            for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {  
           j=out->pattern->index[iptr];  
           /*This can be replaced by bsearch!!*/  
           for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {  
              k=P->pattern->index[jptr];  
              if(k==i) {  
             out->val[iptr]=P->val[jptr];  
             }  
             }  
             }  
             }  
             } else if (block_size==2) {  
                for (i=0;i<out->numRows;++i) {  
                   for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {  
                  j=out->pattern->index[iptr];  
                  /*This can be replaced by bsearch!!*/  
                  for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {  
                     k=P->pattern->index[jptr];  
                     if(k==i) {  
                        out->val[iptr*block_size*block_size]=P->val[jptr*block_size*block_size];  
                        out->val[iptr*block_size*block_size+1]=P->val[jptr*block_size*block_size+2];  
                        out->val[iptr*block_size*block_size+2]=P->val[jptr*block_size*block_size+1];  
                        out->val[iptr*block_size*block_size+3]=P->val[jptr*block_size*block_size+3];  
                        }  
                        }  
                        }  
                        }  
                        } else if (block_size==3) {  
                       for (i=0;i<out->numRows;++i) {  
                          for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {  
                         j=out->pattern->index[iptr];  
                         /*This can be replaced by bsearch!!*/  
                         for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {  
                            k=P->pattern->index[jptr];  
                            if(k==i) {  
                               out->val[iptr*block_size*block_size]=P->val[jptr*block_size*block_size];  
                               out->val[iptr*block_size*block_size+1]=P->val[jptr*block_size*block_size+3];  
                               out->val[iptr*block_size*block_size+2]=P->val[jptr*block_size*block_size+6];  
                               out->val[iptr*block_size*block_size+3]=P->val[jptr*block_size*block_size+1];  
                               out->val[iptr*block_size*block_size+4]=P->val[jptr*block_size*block_size+4];  
                               out->val[iptr*block_size*block_size+5]=P->val[jptr*block_size*block_size+7];  
                               out->val[iptr*block_size*block_size+6]=P->val[jptr*block_size*block_size+2];  
                               out->val[iptr*block_size*block_size+7]=P->val[jptr*block_size*block_size+5];  
                               out->val[iptr*block_size*block_size+8]=P->val[jptr*block_size*block_size+8];  
                               }  
                               }  
                               }  
                               }  
                               }  
                                 
                               /* clean up */  
                               Paso_IndexListArray_free(index_list);  
                               Paso_Pattern_free(outpattern);  
                               return out;  
                               }  
                                 
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    }                            

Legend:
Removed from v.3283  
changed lines
  Added in v.3303

  ViewVC Help
Powered by ViewVC 1.1.26