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    
00054   int iam = In->Comm().MyPID() ;
00055 
00056   int numentries = In->NumGlobalNonzeros();
00057   int NumRowEntries = 0;
00058   double *RowValues = 0;
00059   int *ColIndices = 0;
00060 
00061   int numrows = In->NumGlobalRows();
00062   int numcols = In->NumGlobalCols();
00063 
00064   std::vector <int> Ap( numcols+1 );       // Column i is stored in Aval(Ap[i]..Ap[i+1]-1)
00065   std::vector <int> nextAp( numcols+1 );   // Where to store next value in Column i
00066   std::vector <int> Ai( EPETRA_MAX( numcols, numentries) ) ; //  Row indices
00067   std::vector <double> Aval( EPETRA_MAX( numcols, numentries) ) ; 
00068 
00069   if ( iam == 0 ) { 
00070 
00071     assert( In->NumMyRows() == In->NumGlobalRows() ) ; 
00072     //
00073     //  Count the number of entries in each column
00074     //
00075     std::vector <int>RowsPerCol( numcols ) ; 
00076     for ( int i = 0 ; i < numcols ; i++ ) RowsPerCol[i] = 0 ; 
00077     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00078       assert( In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices ) == 0 ) ;
00079       for ( int j = 0; j < NumRowEntries; j++ ) { 
00080   RowsPerCol[ ColIndices[j] ] ++ ; 
00081       }
00082     }
00083     //
00084     //  Set Ap and nextAp based on RowsPerCol
00085     //
00086     Ap[0] = 0 ; 
00087     for ( int i = 0 ; i < numcols ; i++ ) {
00088       Ap[i+1]= Ap[i] + RowsPerCol[i] ; 
00089       nextAp[i] = Ap[i];
00090     }
00091     //
00092     //  Populate Ai and Aval 
00093     //
00094     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00095       assert( In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices ) == 0 ) ;
00096       for ( int j = 0; j < NumRowEntries; j++ ) { 
00097   Ai[ nextAp[ ColIndices[j] ] ] = MyRow ; 
00098   Aval[ nextAp[ ColIndices[j] ] ] = RowValues[j] ; 
00099   nextAp[ ColIndices[j] ] ++ ; 
00100       }
00101     }
00102 
00103     //
00104     //  Insert values into Out 
00105     //
00106     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00107       int NumInCol = Ap[MyRow+1] -  Ap[MyRow] ;
00108       Out->InsertGlobalValues( MyRow, NumInCol, &Aval[Ap[MyRow]], 
00109          &Ai[Ap[MyRow]] );
00110       assert( Out->IndicesAreGlobal() ) ; 
00111     }
00112   } else {
00113     assert( In->NumMyRows() == 0 ) ; 
00114   }
00115 
00116 
00117   assert( Out->FillComplete()==0 ) ;
00118   return 0 ; 
00119 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines