Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
CrsMatrixTranspose.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29//
30// CrsMatrixTranspose( Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out )
31// fills the matrix "Out" with the transpose of the matrix "In".
32//
33// Notes:
34// Only works on pseudo-distributed matrices, where all rows are stored
35// on process 0.
36// Testing has been limited to square matrices
37// This implementation requires an extra copy of the matrix as an
38// intermediate.
39//
40// One of the bizarre aspects of this code is that it uses the
41// results of calls to ExtractMyRowView() to populate the input
42// to calls to InsertGlobalValues(). This only works because
43// the code is only designed for matrices that all stored entirely on
44// process 0.
45//
46//
47#include "CrsMatrixTranspose.h"
48#include <vector>
49#include "Epetra_Comm.h"
50
51int CrsMatrixTranspose( Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out ) {
52
53#ifndef NDEBUG
54 int ierr = 0;
55#endif
56 int iam = In->Comm().MyPID() ;
57
58 long long numentries = In->NumGlobalNonzeros64();
59 int NumRowEntries = 0;
60 double *RowValues = 0;
61 int *ColIndices = 0;
62
63 long long numrows = In->NumGlobalRows64();
64 long long numcols = In->NumGlobalCols64();
65
66 std::vector <int> Ap( numcols+1 ); // Column i is stored in Aval(Ap[i]..Ap[i+1]-1)
67 std::vector <int> nextAp( numcols+1 ); // Where to store next value in Column i
68#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
69 std::vector <long long> Ai( EPETRA_MAX( numcols, numentries) ) ; // Row indices
70#else
71 std::vector <int> Ai( EPETRA_MAX( numcols, numentries) ) ; // Row indices
72#endif
73 std::vector <double> Aval( EPETRA_MAX( numcols, numentries) ) ;
74
75 if ( iam == 0 ) {
76
77 assert( In->NumMyRows() == In->NumGlobalRows64() ) ;
78 //
79 // Count the number of entries in each column
80 //
81 std::vector <int>RowsPerCol( numcols ) ;
82 for ( int i = 0 ; i < numcols ; i++ ) RowsPerCol[i] = 0 ;
83 for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
84#ifndef NDEBUG
85 ierr =
86#endif
87 In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices );
88 assert( ierr == 0 ) ;
89 for ( int j = 0; j < NumRowEntries; j++ ) {
90 RowsPerCol[ ColIndices[j] ] ++ ;
91 }
92 }
93 //
94 // Set Ap and nextAp based on RowsPerCol
95 //
96 Ap[0] = 0 ;
97 for ( int i = 0 ; i < numcols ; i++ ) {
98 Ap[i+1]= Ap[i] + RowsPerCol[i] ;
99 nextAp[i] = Ap[i];
100 }
101 //
102 // Populate Ai and Aval
103 //
104 for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
105#ifndef NDEBUG
106 ierr =
107#endif
108 In->ExtractMyRowView( MyRow, NumRowEntries, RowValues, ColIndices );
109 assert( ierr == 0 ) ;
110 for ( int j = 0; j < NumRowEntries; j++ ) {
111 Ai[ nextAp[ ColIndices[j] ] ] = MyRow ;
112 Aval[ nextAp[ ColIndices[j] ] ] = RowValues[j] ;
113 nextAp[ ColIndices[j] ] ++ ;
114 }
115 }
116
117 //
118 // Insert values into Out
119 //
120 for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
121 int NumInCol = Ap[MyRow+1] - Ap[MyRow] ;
122#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
123 Out->InsertGlobalValues( MyRow, NumInCol, &Aval[Ap[MyRow]],
124 &Ai[Ap[MyRow]] );
125#endif
126 assert( Out->IndicesAreGlobal() ) ;
127 }
128 } else {
129 assert( In->NumMyRows() == 0 ) ;
130 }
131
132#ifndef NDEBUG
133 ierr =
134#endif
135 Out->FillComplete();
136 assert( ierr==0 ) ;
137 return 0 ;
138}
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)