Amesos Package Browser (Single Doxygen Collection) Development
CrsMatrixTranspose.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 //
00030 //  CrsMatrixTranspose( Epetra_CrsMatrix *In,  Epetra_CrsMatrix *Out ) 
00031 //  fills the matrix "Out" with the transpose of the matrix "In".
00032 //
00033 //  Notes:
00034 //    Only works on pseudo-distributed matrices, where all rows are stored
00035 //    on process 0.
00036 //    Testing has been limited to square matrices
00037 //    This implementation requires an extra copy of the matrix as an 
00038 //    intermediate. 
00039 //
00040 //  One of the bizarre aspects of this code is that it uses the 
00041 //  results of calls to ExtractMyRowView() to populate the input 
00042 //  to calls to InsertGlobalValues().  This only works because
00043 //  the code is only designed for matrices that all stored entirely on 
00044 //  process 0.
00045 //
00046 //
00047 #include "CrsMatrixTranspose.h"
00048 #include <vector>
00049 #include "Epetra_Comm.h"
00050 
00051 int CrsMatrixTranspose( Epetra_CrsMatrix *In,  Epetra_CrsMatrix *Out ) { 
00052 
00053   int ierr = 0;
00054    
00055   int iam = In->Comm().MyPID() ;
00056 
00057   int numentries = In->NumGlobalNonzeros();
00058   int NumRowEntries = 0;
00059   double *RowValues = 0;
00060   int *ColIndices = 0;
00061 
00062   int numrows = In->NumGlobalRows();
00063   int numcols = In->NumGlobalCols();
00064 
00065   std::vector <int> Ap( numcols+1 );       // Column i is stored in Aval(Ap[i]..Ap[i+1]-1)
00066   std::vector <int> nextAp( numcols+1 );   // Where to store next value in Column i
00067   std::vector <int> Ai( EPETRA_MAX( numcols, numentries) ) ; //  Row indices
00068   std::vector <double> Aval( EPETRA_MAX( numcols, numentries) ) ; 
00069 
00070   if ( iam == 0 ) { 
00071 
00072     assert( In->NumMyRows() == In->NumGlobalRows() ) ; 
00073     //
00074     //  Count the number of entries in each column
00075     //
00076     std::vector <int>RowsPerCol( numcols ) ; 
00077     for ( int i = 0 ; i < numcols ; i++ ) RowsPerCol[i] = 0 ; 
00078     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00079       ierr = In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices );
00080       assert( ierr == 0 ) ;
00081       for ( int j = 0; j < NumRowEntries; j++ ) { 
00082   RowsPerCol[ ColIndices[j] ] ++ ; 
00083       }
00084     }
00085     //
00086     //  Set Ap and nextAp based on RowsPerCol
00087     //
00088     Ap[0] = 0 ; 
00089     for ( int i = 0 ; i < numcols ; i++ ) {
00090       Ap[i+1]= Ap[i] + RowsPerCol[i] ; 
00091       nextAp[i] = Ap[i];
00092     }
00093     //
00094     //  Populate Ai and Aval 
00095     //
00096     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00097       ierr = In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices );
00098       assert( ierr == 0 ) ;
00099       for ( int j = 0; j < NumRowEntries; j++ ) { 
00100   Ai[ nextAp[ ColIndices[j] ] ] = MyRow ; 
00101   Aval[ nextAp[ ColIndices[j] ] ] = RowValues[j] ; 
00102   nextAp[ ColIndices[j] ] ++ ; 
00103       }
00104     }
00105 
00106     //
00107     //  Insert values into Out 
00108     //
00109     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00110       int NumInCol = Ap[MyRow+1] -  Ap[MyRow] ;
00111       Out->InsertGlobalValues( MyRow, NumInCol, &Aval[Ap[MyRow]], 
00112          &Ai[Ap[MyRow]] );
00113       assert( Out->IndicesAreGlobal() ) ; 
00114     }
00115   } else {
00116     assert( In->NumMyRows() == 0 ) ; 
00117   }
00118 
00119   ierr = Out->FillComplete();
00120   assert( ierr==0 ) ;
00121   return 0 ; 
00122 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines