1 /* _______________________________________________________________________
2
3 PECOS: Parallel Environment for Creation Of Stochastics
4 Copyright (c) 2011, Sandia National Laboratories.
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Pecos directory.
7 _______________________________________________________________________ */
8
9 #include "sparse_matrices.hpp"
10
11 namespace Pecos {
12 namespace util {
13
get_row(int row_num,RealMatrix & row) const14 void BlockDiagonalMatrix::get_row( int row_num, RealMatrix &row ) const {
15 int num_total_rows = 0;
16 int num_total_rows_prev = 0;
17 int block_row = 0;
18 int i=0;
19 for (i=0; i<numBlocks_; i++ ){
20 num_total_rows_prev = num_total_rows;
21 num_total_rows += blocks_[i].numRows();
22 if ( row_num < num_total_rows ){
23 block_row = row_num - num_total_rows_prev;
24 break;
25 }
26 }
27 if ( ( row.numRows() !=1 ) || (row.numCols() != blocks_[i].numCols() ) )
28 row.shapeUninitialized( 1, blocks_[i].numCols() );
29 for ( int j=0; j<blocks_[i].numCols(); j++)
30 row(0,j) = blocks_[i](block_row,j);
31 }
32
pre_multiply(const RealMatrix & matrix,RealMatrix & result,Teuchos::ETransp block_trans) const33 void BlockDiagonalMatrix::pre_multiply( const RealMatrix &matrix,
34 RealMatrix &result,
35 Teuchos::ETransp block_trans ) const {
36 int block_num_cols = num_cols();
37 int result_num_rows = num_rows();
38 if ( block_trans==Teuchos::TRANS ){
39 block_num_cols = num_rows();
40 result_num_rows = num_cols();
41 }
42
43 if ( block_num_cols != matrix.numRows() ){
44 std::string msg = "BlockDiagonalMatrix::pre_multiply() ";
45 msg += "Matrices sizes are inconsistent\n";
46 throw( std::runtime_error( msg ) );
47 }
48
49 result.shapeUninitialized( result_num_rows, matrix.numCols() );
50 int sub_matrix_start_row = 0;
51 int sub_result_start_row = 0;
52 for (int i=0; i<numBlocks_; i++ ){
53 int num_block_rows = blocks_[i].numRows();
54 int num_block_cols = blocks_[i].numCols();
55 if ( block_trans == Teuchos::TRANS ){
56 num_block_rows = blocks_[i].numCols();
57 num_block_cols = blocks_[i].numRows();
58 }
59 int num_submatrix_rows = num_block_cols;
60 RealMatrix sub_matrix( Teuchos::View, matrix,
61 num_submatrix_rows, matrix.numCols(),
62 sub_matrix_start_row, 0 );
63 RealMatrix sub_result( Teuchos::View, result,
64 num_block_rows, matrix.numCols(),
65 sub_result_start_row, 0 );
66 sub_result.multiply( block_trans, Teuchos::NO_TRANS,
67 1.0, blocks_[i], sub_matrix, 0.0 );
68 sub_matrix_start_row += num_submatrix_rows;
69 sub_result_start_row += num_block_rows;
70 }
71 }
72
73
post_multiply(const RealMatrix & matrix,RealMatrix & result,Teuchos::ETransp block_trans,int subblock_num_rows) const74 void BlockDiagonalMatrix::post_multiply( const RealMatrix &matrix,
75 RealMatrix &result,
76 Teuchos::ETransp block_trans,
77 int subblock_num_rows ) const {
78 int block_num_rows = num_rows();
79 int result_num_cols = num_cols();
80 if ( block_trans==Teuchos::TRANS ){
81 block_num_rows = num_cols();
82 result_num_cols = num_rows();
83 }
84
85 if ( subblock_num_rows < 0 )
86 subblock_num_rows = block_num_rows;
87
88 if ( subblock_num_rows != matrix.numCols() ){
89 std::string msg = "BlockDiagonalMatrix::post_multiply() ";
90 msg += "Matrices sizes are inconsistent\n";
91 throw( std::runtime_error( msg ) );
92 }
93
94 if ( subblock_num_rows > block_num_rows ){
95 std::string msg = "BlockDiagonalMatrix::post_multiply() ";
96 msg += "The number of subset rows (subblock_num_rows) is to large\n";
97 throw( std::runtime_error( msg ) );
98 }
99
100
101 result.shape( matrix.numRows(), result_num_cols );
102 int sub_matrix_start_col = 0;
103 int sub_result_start_col = 0;
104
105 bool done = false;
106 int num_total_rows = 0;
107 int num_total_rows_prev = 0;
108 for (int i=0; i<numBlocks_; i++ ){
109 int num_block_rows = blocks_[i].numRows();
110 int num_block_cols = blocks_[i].numCols();
111 if ( block_trans == Teuchos::TRANS ){
112 num_block_rows = blocks_[i].numCols();
113 num_block_cols = blocks_[i].numRows();
114 }
115
116 num_total_rows_prev = num_total_rows;
117 num_total_rows += num_block_rows;
118 int num_sub_block_rows = num_block_rows;
119 int num_sub_block_cols = num_block_cols;
120 if ( num_total_rows > subblock_num_rows ){
121 num_sub_block_rows = subblock_num_rows - num_total_rows_prev;
122 done = true;
123 }
124 int num_no_trans_sub_block_rows = num_sub_block_rows;
125 int num_no_trans_sub_block_cols = num_sub_block_cols;
126 if ( block_trans == Teuchos::TRANS ){
127 num_no_trans_sub_block_rows = num_sub_block_cols;
128 num_no_trans_sub_block_cols = num_sub_block_rows;
129 }
130 RealMatrix sub_block( Teuchos::View, blocks_[i], num_no_trans_sub_block_rows,
131 num_no_trans_sub_block_cols, 0, 0 );
132
133 int num_submatrix_cols = num_sub_block_rows;
134 RealMatrix sub_matrix( Teuchos::View, matrix,
135 matrix.numRows(), num_submatrix_cols,
136 0, sub_matrix_start_col );
137 RealMatrix sub_result( Teuchos::View, result,
138 matrix.numRows(), num_sub_block_cols,
139 0, sub_result_start_col );
140 sub_result.multiply( Teuchos::NO_TRANS, block_trans,
141 //1.0, sub_matrix, blocks_[i], 0.0 );
142 1.0, sub_matrix, sub_block, 0.0 );
143
144 sub_matrix_start_col += num_submatrix_cols;
145 sub_result_start_col += num_sub_block_cols;
146
147 if ( done ) break;
148 }
149 }
150
post_multiply_block(int block_num,const RealMatrix & matrix,RealMatrix & result) const151 void BlockDiagonalMatrix::post_multiply_block( int block_num,
152 const RealMatrix &matrix,
153 RealMatrix & result ) const
154 {
155 if ( block_num >= num_blocks() ){
156 std::string msg = "BlockDiagonalMatrix::post_multiply_block() ";
157 msg += "block num exceeds the number of blocks\n";
158 throw( std::runtime_error( msg ) );
159 }
160
161 int block_num_rows = blocks_[block_num].numRows();
162 int block_num_cols = blocks_[block_num].numCols();
163
164 int matrix_num_rows = matrix.numRows();
165 int matrix_num_cols = matrix.numCols();
166
167 if ( block_num_rows != matrix_num_cols ){
168 std::string msg = "BlockDiagonalMatrix::post_multiply_block() ";
169 msg += "Matrices sizes are inconsistent\n";
170 throw( std::runtime_error( msg ) );
171 }
172
173 result.shapeUninitialized( matrix_num_rows, block_num_cols );
174 result.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS,
175 1.0, matrix, blocks_[block_num], 0.0 );
176 };
177
178 } // namespace util
179 } // namespace Pecos
180