This package works
This commit is contained in:
vasiloglou
Executable
+40
@@ -0,0 +1,40 @@
|
||||
1 1 0.31106
|
||||
13 1 0.47969
|
||||
2 2 0.734
|
||||
3 3 0.4109
|
||||
4 3 0.6412
|
||||
1 4 0.16853
|
||||
3 4 0.39979
|
||||
11 4 0.87535
|
||||
19 4 0.9667
|
||||
20 4 0.31775
|
||||
14 5 0.61591
|
||||
1 6 0.89665
|
||||
15 6 0.61663
|
||||
11 7 0.83526
|
||||
1 8 0.32272
|
||||
6 8 0.69778
|
||||
18 8 0.51015
|
||||
3 9 0.50552
|
||||
6 9 0.46189
|
||||
18 9 0.71396
|
||||
13 10 0.56082
|
||||
14 11 0.6619
|
||||
20 11 0.5877
|
||||
3 12 0.16931
|
||||
6 12 0.082613
|
||||
7 12 0.82072
|
||||
5 13 0.83685
|
||||
15 13 0.68514
|
||||
8 14 0.19302
|
||||
18 14 0.51521
|
||||
12 15 0.88071
|
||||
5 16 0.80346
|
||||
8 16 0.44535
|
||||
10 16 0.30874
|
||||
3 17 0.52475
|
||||
4 18 0.016197
|
||||
11 18 0.3331
|
||||
19 18 0.82212
|
||||
18 19 0.60587
|
||||
9 20 0.012958
|
||||
Executable
+6
@@ -0,0 +1,6 @@
|
||||
3 3 0.37376
|
||||
19 4 0.04147
|
||||
14 5 0.1847
|
||||
1 6 0.20795
|
||||
20 11 0.43719
|
||||
18 14 0.13314
|
||||
Executable
+73
@@ -0,0 +1,73 @@
|
||||
1 1 0.31106
|
||||
13 1 0.47969
|
||||
16 1 -0.11207
|
||||
2 2 0.734
|
||||
3 3 -0.4987
|
||||
4 3 0.6412
|
||||
12 3 -0.89951
|
||||
16 3 -0.29156
|
||||
1 4 0.16853
|
||||
3 4 0.39979
|
||||
10 4 -0.67427
|
||||
11 4 0.87535
|
||||
19 4 0.9238
|
||||
20 4 0.31775
|
||||
10 5 -0.9271
|
||||
14 5 0.31602
|
||||
19 5 -0.0058848
|
||||
1 6 0.66473
|
||||
8 6 -0.199
|
||||
15 6 0.61663
|
||||
17 6 -0.94423
|
||||
11 7 0.83526
|
||||
12 7 -0.69276
|
||||
16 7 -0.097447
|
||||
1 8 0.32272
|
||||
6 8 0.69778
|
||||
18 8 0.51015
|
||||
3 9 0.50552
|
||||
6 9 0.46189
|
||||
10 9 -0.34382
|
||||
16 9 -0.39745
|
||||
18 9 0.71396
|
||||
1 10 -0.47866
|
||||
5 10 -0.76755
|
||||
13 10 0.56082
|
||||
1 11 -0.52652
|
||||
14 11 0.6619
|
||||
19 11 -0.57442
|
||||
20 11 -0.15621
|
||||
3 12 0.16931
|
||||
6 12 0.082613
|
||||
7 12 0.82072
|
||||
10 12 -0.59449
|
||||
13 12 -0.60971
|
||||
16 12 -0.33331
|
||||
5 13 0.83685
|
||||
6 13 -0.94734
|
||||
15 13 0.68514
|
||||
3 14 -0.9222
|
||||
6 14 -0.81331
|
||||
7 14 -0.92383
|
||||
8 14 0.19302
|
||||
18 14 0.25678
|
||||
1 15 -0.79272
|
||||
12 15 0.88071
|
||||
2 16 -0.19301
|
||||
5 16 0.80346
|
||||
8 16 0.44535
|
||||
10 16 0.30874
|
||||
11 16 -0.0033741
|
||||
14 16 -0.85604
|
||||
3 17 0.52475
|
||||
4 18 0.016197
|
||||
11 18 0.3331
|
||||
12 18 -0.43965
|
||||
19 18 0.82212
|
||||
4 19 -0.013266
|
||||
11 19 -0.98201
|
||||
17 19 -0.83856
|
||||
18 19 0.60587
|
||||
9 20 0.012958
|
||||
10 20 -0.61549
|
||||
12 20 -0.70102
|
||||
Executable
+73
@@ -0,0 +1,73 @@
|
||||
1 1 0.31106
|
||||
13 1 0.47969
|
||||
16 1 0.11207
|
||||
2 2 0.734
|
||||
3 3 1.3205
|
||||
4 3 0.6412
|
||||
12 3 0.89951
|
||||
16 3 0.29156
|
||||
1 4 0.16853
|
||||
3 4 0.39979
|
||||
10 4 0.67427
|
||||
11 4 0.87535
|
||||
19 4 1.0096
|
||||
20 4 0.31775
|
||||
10 5 0.9271
|
||||
14 5 0.91579
|
||||
19 5 0.0058848
|
||||
1 6 1.1286
|
||||
8 6 0.199
|
||||
15 6 0.61663
|
||||
17 6 0.94423
|
||||
11 7 0.83526
|
||||
12 7 0.69276
|
||||
16 7 0.097447
|
||||
1 8 0.32272
|
||||
6 8 0.69778
|
||||
18 8 0.51015
|
||||
3 9 0.50552
|
||||
6 9 0.46189
|
||||
10 9 0.34382
|
||||
16 9 0.39745
|
||||
18 9 0.71396
|
||||
1 10 0.47866
|
||||
5 10 0.76755
|
||||
13 10 0.56082
|
||||
1 11 0.52652
|
||||
14 11 0.6619
|
||||
19 11 0.57442
|
||||
20 11 1.3316
|
||||
3 12 0.16931
|
||||
6 12 0.082613
|
||||
7 12 0.82072
|
||||
10 12 0.59449
|
||||
13 12 0.60971
|
||||
16 12 0.33331
|
||||
5 13 0.83685
|
||||
6 13 0.94734
|
||||
15 13 0.68514
|
||||
3 14 0.9222
|
||||
6 14 0.81331
|
||||
7 14 0.92383
|
||||
8 14 0.19302
|
||||
18 14 0.77364
|
||||
1 15 0.79272
|
||||
12 15 0.88071
|
||||
2 16 0.19301
|
||||
5 16 0.80346
|
||||
8 16 0.44535
|
||||
10 16 0.30874
|
||||
11 16 0.0033741
|
||||
14 16 0.85604
|
||||
3 17 0.52475
|
||||
4 18 0.016197
|
||||
11 18 0.3331
|
||||
12 18 0.43965
|
||||
19 18 0.82212
|
||||
4 19 0.013266
|
||||
11 19 0.98201
|
||||
17 19 0.83856
|
||||
18 19 0.60587
|
||||
9 20 0.012958
|
||||
10 20 0.61549
|
||||
12 20 0.70102
|
||||
Executable
+70
@@ -0,0 +1,70 @@
|
||||
5 1 0.090046
|
||||
8 1 0.049912
|
||||
10 1 0.034601
|
||||
3 3 0.52605
|
||||
4 3 0.58324
|
||||
5 3 0.23426
|
||||
6 3 0.074311
|
||||
7 3 0.73825
|
||||
8 3 0.12985
|
||||
10 3 0.090017
|
||||
13 4 0.37814
|
||||
18 4 0.025991
|
||||
8 5 0.057884
|
||||
13 5 0.51993
|
||||
18 5 0.15807
|
||||
1 6 0.13636
|
||||
3 6 0.49548
|
||||
6 6 0.13886
|
||||
13 6 0.11125
|
||||
18 6 0.10152
|
||||
3 7 0.11729
|
||||
5 7 0.078295
|
||||
6 7 0.05723
|
||||
7 7 0.56856
|
||||
8 7 0.043398
|
||||
10 7 0.030086
|
||||
5 9 0.31933
|
||||
8 9 0.17701
|
||||
10 9 0.12271
|
||||
13 9 0.19282
|
||||
1 10 0.14889
|
||||
13 10 0.22961
|
||||
14 10 0.47274
|
||||
1 11 0.16378
|
||||
9 11 0.0096394
|
||||
13 11 0.25257
|
||||
18 11 0.34802
|
||||
5 12 0.77805
|
||||
8 12 0.14844
|
||||
10 12 0.10291
|
||||
13 12 0.3334
|
||||
15 12 0.41774
|
||||
1 13 0.84943
|
||||
15 13 0.58416
|
||||
1 14 0.72925
|
||||
3 14 0.37893
|
||||
4 14 0.5955
|
||||
11 14 0.85772
|
||||
15 14 0.50151
|
||||
19 14 0.21246
|
||||
1 15 0.24658
|
||||
13 15 0.38026
|
||||
2 16 0.14167
|
||||
8 16 0.16523
|
||||
14 16 0.0022333
|
||||
18 16 0.44104
|
||||
20 16 0.0019829
|
||||
3 18 0.074436
|
||||
6 18 0.036321
|
||||
7 18 0.36083
|
||||
1 19 0.0022358
|
||||
3 19 0.44533
|
||||
11 19 0.011613
|
||||
14 19 0.64999
|
||||
19 19 0.012825
|
||||
20 19 0.58134
|
||||
3 20 0.11869
|
||||
6 20 0.057913
|
||||
7 20 0.57534
|
||||
13 20 0.34518
|
||||
Executable
+39
@@ -0,0 +1,39 @@
|
||||
16 1 0.11207
|
||||
3 3 0.9096
|
||||
12 3 0.89951
|
||||
16 3 0.29156
|
||||
10 4 0.67427
|
||||
19 4 0.042898
|
||||
10 5 0.9271
|
||||
14 5 0.29989
|
||||
19 5 0.0058848
|
||||
1 6 0.23192
|
||||
8 6 0.199
|
||||
17 6 0.94423
|
||||
12 7 0.69276
|
||||
16 7 0.097447
|
||||
10 9 0.34382
|
||||
16 9 0.39745
|
||||
1 10 0.47866
|
||||
5 10 0.76755
|
||||
1 11 0.52652
|
||||
19 11 0.57442
|
||||
20 11 0.7439
|
||||
10 12 0.59449
|
||||
13 12 0.60971
|
||||
16 12 0.33331
|
||||
6 13 0.94734
|
||||
3 14 0.9222
|
||||
6 14 0.81331
|
||||
7 14 0.92383
|
||||
18 14 0.25843
|
||||
1 15 0.79272
|
||||
2 16 0.19301
|
||||
11 16 0.0033741
|
||||
14 16 0.85604
|
||||
12 18 0.43965
|
||||
4 19 0.013266
|
||||
11 19 0.98201
|
||||
17 19 0.83856
|
||||
10 20 0.61549
|
||||
12 20 0.70102
|
||||
@@ -0,0 +1,12 @@
|
||||
librule(name="sparse",
|
||||
headers=["sparse_matrix.h"],
|
||||
# ]
|
||||
deplibs=["trilinos:trilinos", "base:base"]
|
||||
);
|
||||
|
||||
binrule(name="mtest",
|
||||
sources=["sparse_matrix_test.cc"],
|
||||
cflags=" -fexceptions",
|
||||
deplibs=[":sparse"]
|
||||
);
|
||||
|
||||
@@ -0,0 +1,550 @@
|
||||
/*
|
||||
* =====================================================================================
|
||||
*
|
||||
* Filename: sparse_matrix.h
|
||||
*
|
||||
* Description:
|
||||
*
|
||||
* Version: 1.0
|
||||
* Created: 12/01/2007 04:12:00 PM EST
|
||||
* Revision: none
|
||||
* Compiler: gcc
|
||||
*
|
||||
* Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org
|
||||
* Company: Georgia Tech Fastlab-ESP Lab
|
||||
*
|
||||
* =====================================================================================
|
||||
*/
|
||||
|
||||
#ifndef SPARSE_MATRIX_H_
|
||||
#define SPARSE_MATRIX_H_
|
||||
#ifndef HAVE_CONFIG_H
|
||||
#define HAVE_CONFIG_H
|
||||
#endif
|
||||
#ifndef USE_TRILINOS
|
||||
#define USE_TRILINOS
|
||||
#endif
|
||||
#include <stdio.h>
|
||||
#include <errno.h>
|
||||
#include <string>
|
||||
#include <map>
|
||||
#include <vector>
|
||||
#include <sstream>
|
||||
#include <algorithm>
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "la/matrix.h"
|
||||
// you need this because trillinos redifines it. It's ok
|
||||
// if you don't have it, but you will get an annoying warning
|
||||
#ifdef F77_FUNC
|
||||
#undef F77_FUNC
|
||||
#endif
|
||||
#include "trilinos/include/Epetra_CrsMatrix.h"
|
||||
#include "trilinos/include/Epetra_SerialComm.h"
|
||||
#include "trilinos/include/Epetra_Map.h"
|
||||
#include "trilinos/include/Epetra_Vector.h"
|
||||
#include "trilinos/include/Epetra_MultiVector.h"
|
||||
#include "trilinos/include/AnasaziBasicEigenproblem.hpp"
|
||||
#include "trilinos/include/AnasaziEpetraAdapter.hpp"
|
||||
#include "trilinos/include/AnasaziBlockKrylovSchurSolMgr.hpp"
|
||||
#include "trilinos/include/AztecOO.h"
|
||||
|
||||
/* class SparseMatrix created by Nick
|
||||
* This is a sparse matrix wrapper for trilinos Epetra_CrsMatrix
|
||||
* It is much simpler than Epetra_CrsMatrix. At this time
|
||||
* it supports eigenvalues (Krylov method) and linear system solution
|
||||
* I have added matrix addition/subtraction multiplication
|
||||
* I am also trying to add the submatrices
|
||||
* Note: There is a restriction on these matrices, the number of rows is
|
||||
* always greater or equal to the number of columns. The number of rows
|
||||
* is also called dimension. We pose this restriction because trilinos supports
|
||||
* square matrices only. In sparse matrices though this is not the problem since
|
||||
* an mxn matrix where m>n can is equivalent to an mxm matrix where all the
|
||||
* elements with n<j<m are zero
|
||||
*/
|
||||
|
||||
class Sparsem;
|
||||
class SparseMatrix {
|
||||
public:
|
||||
friend class Sparsem;
|
||||
// Some typedefs for oft-used data types
|
||||
typedef Epetra_MultiVector MV;
|
||||
typedef Epetra_Operator OP;
|
||||
typedef Anasazi::MultiVecTraits<double, Epetra_MultiVector> MVT;
|
||||
|
||||
SparseMatrix() ;
|
||||
// Constructor
|
||||
// num_of_rows: number of rows
|
||||
// num_of_cols: number of columns
|
||||
// nnz_per_row: an estimate of the non zero elements per row
|
||||
// This doesn't need to be accurate. If you need
|
||||
// more it will automatically resize. Try to be as accurate
|
||||
// as you can because resizing costs. It is better if your
|
||||
// estimete if greater than the true non zero elements. So
|
||||
// it is better to overestimate than underestimate
|
||||
SparseMatrix(const index_t num_of_rows,
|
||||
const index_t num_of_cols,
|
||||
const index_t nnz_per_row);
|
||||
// Copy constructor
|
||||
SparseMatrix(const SparseMatrix &other);
|
||||
SparseMatrix(std::string textfile) {
|
||||
Init(textfile);
|
||||
}
|
||||
~SparseMatrix() {
|
||||
Destruct();
|
||||
}
|
||||
void Destruct();
|
||||
// Use this initializer like the Constructor
|
||||
void Init(const index_t num_of_rows,
|
||||
const index_t num_of_columns,
|
||||
const index_t nnz_per_row);
|
||||
// This Initializer is like the previous one with the main difference that
|
||||
// for every row we give a seperate estimate for the non-zero elements.
|
||||
void Init(index_t num_of_rows, index_t num_of_columns, index_t *nnz_per_row);
|
||||
// This Initializer fills the sparse matrix with data.
|
||||
// row_indices: row indices for non-zero elements
|
||||
// col_indices: column indices for non-zero elements
|
||||
// values : values of non-zeros elements
|
||||
// If the dimension (number of rows)and the expected (nnz elements per row)
|
||||
// are set to a negative value, the function will automatically detect it
|
||||
void Init(const std::vector<index_t> &row_indices,
|
||||
const std::vector<index_t> &col_indices,
|
||||
const Vector &values,
|
||||
index_t nnz_per_row,
|
||||
index_t dimension);
|
||||
// The same as above but we use STL vector for values
|
||||
void Init(const std::vector<index_t> &row_indices,
|
||||
const std::vector<index_t> &col_indices,
|
||||
const std::vector<double> &values,
|
||||
index_t nnz_per_row,
|
||||
index_t dimension);
|
||||
|
||||
// Initialize from a text file in the following format
|
||||
// row column value \n
|
||||
void Init(std::string textfile);
|
||||
// Initialize the diagonal
|
||||
void InitDiagonal(const Vector &vec);
|
||||
// Initialize the diagonal with a constant
|
||||
void InitDiagonal(const double value);
|
||||
// It is recomended that you load the matrix row-wise, Before
|
||||
// you do that call StartLoadingRows()
|
||||
void StartLoadingRows();
|
||||
// All these functions load Rows, with the data in different format
|
||||
void LoadRow(index_t row, std::vector<index_t> &columns, Vector &values);
|
||||
void LoadRow(index_t row, index_t *columns, Vector &values);
|
||||
void LoadRow(index_t row, index_t num, index_t *columns, double *values);
|
||||
void LoadRow(index_t row, std::vector<index_t> &columns, std::vector<double> &values);
|
||||
// When you are done call this it does some optimization in the storage, no
|
||||
// further asignment
|
||||
void EndLoading();
|
||||
// It makes the matrix symmetric. It scans the rows of the matrix and for every (i,j)
|
||||
// element (j,i) equal to (j,i)
|
||||
void MakeSymmetric();
|
||||
// if you know that the matrix is symmetric set the flag
|
||||
void set_symmetric(bool val) {
|
||||
issymmetric_ = val;
|
||||
}
|
||||
// Not implemented yet
|
||||
void SetDiagonal(const Vector &vector);
|
||||
// Copy function, used also by copy constructor
|
||||
void Copy(const SparseMatrix &other);
|
||||
// Not implemented yet
|
||||
void Alias(const Matrix& other);
|
||||
// Not Implemented yet
|
||||
void SwapValues(SparseMatrix* other);
|
||||
// Access values, It will fail if EndLoading() has been called
|
||||
double get(index_t r, index_t c) const;
|
||||
// Set Values
|
||||
void set(index_t r, index_t c, double v);
|
||||
// For debug purposes you can call it to print the matrix
|
||||
std::string Print() {
|
||||
std::ostringstream s1;
|
||||
matrix_->Print(s1);
|
||||
return s1.str();
|
||||
}
|
||||
// Get the number of rows
|
||||
index_t get_num_of_rows() {
|
||||
return num_of_rows_;
|
||||
}
|
||||
// Get the number of columns
|
||||
index_t get_num_of_columns() {
|
||||
return num_of_columns_;
|
||||
}
|
||||
// Dimension should be equal to the number of rows
|
||||
index_t get_dimension() {
|
||||
return dimension_;
|
||||
}
|
||||
// The number of non zero elements
|
||||
index_t get_nnz() {
|
||||
return matrix_->NumGlobalNonzeros();
|
||||
}
|
||||
// Computes the eignvalues with the Krylov Method
|
||||
void Eig(index_t num_of_eigvalues, // number of eigenvalues to compute
|
||||
std::string eigtype, // Choose which eigenvalues to compute
|
||||
// Choices are:
|
||||
// LM - target the largest magnitude
|
||||
// SM - target the smallest magnitude
|
||||
// LR - target the largest real
|
||||
// SR - target the smallest real
|
||||
// LI - target the largest imaginary
|
||||
// SI - target the smallest imaginary
|
||||
Matrix *eigvectors, // The eigenvectors computed must not be initialized
|
||||
std::vector<double> *real_eigvalues, // real part of the eigenvalues
|
||||
// must be initialized, but should not
|
||||
// allocate space. The eigenvalues
|
||||
// returned might actually be less
|
||||
// than the ones requested
|
||||
// for example when the matrix has
|
||||
// rank n< eigenvalues requested
|
||||
std::vector<double> *imag_eigvalues // imaginary part of the eigenvalues
|
||||
// must be initialized. If the
|
||||
// problem is symmetric there is
|
||||
// no need to initialize.
|
||||
// The same as real_eigvalues
|
||||
// hold for the space allocated
|
||||
// in the non-symmetric case
|
||||
);
|
||||
// Linear System solution, Call Endloading First.
|
||||
void LinSolve(Vector &b, // must be initialized (space allocated)
|
||||
Vector *x, // must be initialized (space allocated)
|
||||
double tolerance,
|
||||
index_t iterations
|
||||
) {
|
||||
if (matrix_->Filled()==false && matrix_->StorageOptimized()==false){
|
||||
FATAL("You should call EndLoading() first\n");
|
||||
}
|
||||
Epetra_Vector tempb(View, *map_, b.ptr());
|
||||
Epetra_Vector tempx(View, *map_, x->ptr());
|
||||
// create linear problem
|
||||
Epetra_LinearProblem problem(matrix_.get(), &tempx, &tempb);
|
||||
// create the AztecOO instance
|
||||
AztecOO solver(problem);
|
||||
solver.SetAztecOption( AZ_precond, AZ_Jacobi);
|
||||
solver.Iterate(iterations, tolerance);
|
||||
NONFATAL("Solver performed %i iterations, true residual %lg",
|
||||
solver.NumIters(), solver.TrueResidual());
|
||||
}
|
||||
|
||||
// Use this for the general case
|
||||
void LinSolve(Vector &b, Vector *x) {
|
||||
LinSolve(b, x, 1E-9, 1000);
|
||||
}
|
||||
|
||||
// scales the matrix with a scalar;
|
||||
void Scale(double scalar) {
|
||||
matrix_->Scale(scalar);
|
||||
}
|
||||
// The matrix will be scaled such that A(i,j) = x(j)*A(i,j)
|
||||
// where i denotes the global row number of A and j denotes the column number
|
||||
void ColumnScale(const Vector &vec) {
|
||||
Epetra_Vector temp(View, *map_, (double*)vec.ptr());
|
||||
matrix_->RightScale(temp);
|
||||
}
|
||||
// The matrix will be scaled such that A(i,j) = x(i)*A(i,j)
|
||||
// where i denotes the row number of A and j denotes the column number of A.
|
||||
void RowScale(const Vector &vec) {
|
||||
Epetra_Vector temp(View, *map_, (double *)vec.ptr());
|
||||
matrix_->LeftScale(temp);
|
||||
}
|
||||
// computes the L1 norm
|
||||
double L1Norm() {
|
||||
return matrix_->NormOne();
|
||||
}
|
||||
// L infinity norm
|
||||
double LInfNorm() {
|
||||
return matrix_->NormInf();
|
||||
}
|
||||
// Computes the inverse of the sum of absolute values of the rows
|
||||
// of the matrix
|
||||
void InvRowsSums(Vector *result) {
|
||||
Epetra_Vector temp(View, *map_, result->ptr());
|
||||
matrix_->InvRowSums(temp);
|
||||
}
|
||||
// Computes the inv of max of absolute values of the rows of the matrix,
|
||||
void InvRowMaxs(Vector *result) {
|
||||
Epetra_Vector temp(View, *map_, result->ptr());
|
||||
matrix_->InvRowMaxs(temp);
|
||||
}
|
||||
// Computes the inverse of the sum of absolute values of the columns of the
|
||||
// matrix
|
||||
void InvColSums(Vector *result) {
|
||||
Epetra_Vector temp(View, *map_, result->ptr());
|
||||
matrix_->InvColSums(temp);
|
||||
}
|
||||
// Computes the inv of max of absolute values of the columns of the matrix,
|
||||
void InvColMaxs(Vector *result) {
|
||||
Epetra_Vector temp(View, *map_, result->ptr());
|
||||
matrix_->InvColMaxs(temp);
|
||||
}
|
||||
|
||||
private:
|
||||
index_t dimension_;
|
||||
index_t num_of_rows_;
|
||||
index_t num_of_columns_;
|
||||
Epetra_SerialComm comm_;
|
||||
bool issymmetric_;
|
||||
Epetra_Map *map_;
|
||||
Teuchos::RCP<Epetra_CrsMatrix> matrix_;
|
||||
index_t *my_global_elements_;
|
||||
void Load(const std::vector<index_t> &rows,
|
||||
const std::vector<index_t> &columns,
|
||||
const Vector &values);
|
||||
void AllRowsLoad(Vector &rows, Vector &columns);
|
||||
|
||||
};
|
||||
|
||||
class Sparsem {
|
||||
public:
|
||||
static inline void Add(const SparseMatrix &a,
|
||||
const SparseMatrix &b,
|
||||
SparseMatrix *result) {
|
||||
DEBUG_ASSERT(a.num_of_rows_==b.num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_columns_==b.num_of_columns_);
|
||||
DEBUG_ASSERT(a.num_of_rows_==result->num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_columns_==result->num_of_columns_);
|
||||
result->StartLoadingRows();
|
||||
for(index_t r=0; r<a.num_of_rows_; r++) {
|
||||
index_t num1, num2;
|
||||
double *values1, *values2;
|
||||
index_t *indices1, *indices2;
|
||||
a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1);
|
||||
b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2);
|
||||
std::vector<double> values3;
|
||||
std::vector<index_t> indices3;
|
||||
index_t i=0;
|
||||
index_t j=0;
|
||||
while (likely(i<num1 && j<num2)) {
|
||||
while (indices1[i] < indices2[j]) {
|
||||
values3.push_back(values1[i]);
|
||||
indices3.push_back(indices1[i]);
|
||||
i++;
|
||||
if unlikely((i>=num1)) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
if ( likely(i<num1) && indices1[i] == indices2[j]) {
|
||||
values3.push_back(values1[i] + values2[j]);
|
||||
indices3.push_back(indices1[i]);
|
||||
} else {
|
||||
values3.push_back(values2[j]);
|
||||
indices3.push_back(indices2[j]);
|
||||
}
|
||||
j++;
|
||||
}
|
||||
if (i<num1) {
|
||||
values3.insert(values3.end(), values1+i, values1+num1);
|
||||
indices3.insert(indices3.end(), indices1+i, indices1+num1);
|
||||
}
|
||||
if (j<num2) {
|
||||
values3.insert(values3.end(), values2+j, values2+num2);
|
||||
indices3.insert(indices3.end(), indices2+j, indices2+num2);
|
||||
}
|
||||
result->LoadRow(r, indices3, values3);
|
||||
}
|
||||
}
|
||||
static inline void Subtract(const SparseMatrix &a,
|
||||
const SparseMatrix &b,
|
||||
SparseMatrix *result) {
|
||||
DEBUG_ASSERT(a.num_of_rows_==b.num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_columns_==b.num_of_columns_);
|
||||
DEBUG_ASSERT(a.num_of_rows_==result->num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_columns_==result->num_of_columns_);
|
||||
// If you try assigning the results to an already initialized matrix
|
||||
// you might get unexpected results. The following assertions
|
||||
// prevent you partially from that
|
||||
DEBUG_ASSERT(&a != result);
|
||||
DEBUG_ASSERT(&b != result);
|
||||
result->StartLoadingRows();
|
||||
for(index_t r=0; r<a.num_of_rows_; r++) {
|
||||
index_t num1, num2;
|
||||
double *values1, *values2;
|
||||
index_t *indices1, *indices2;
|
||||
a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1);
|
||||
b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2);
|
||||
std::vector<double> values3;
|
||||
std::vector<index_t> indices3;
|
||||
index_t i=0;
|
||||
index_t j=0;
|
||||
while (likely(i<num1 && j<num2)) {
|
||||
while (indices1[i] < indices2[j]) {
|
||||
values3.push_back(values1[i]);
|
||||
indices3.push_back(indices1[i]);
|
||||
i++;
|
||||
if unlikely((i>=num1)) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (likely(i<num1) && indices1[i] == indices2[j]) {
|
||||
double diff=values1[i] - values2[j];
|
||||
if (diff!=0) {
|
||||
values3.push_back(diff);
|
||||
indices3.push_back(indices1[i]);
|
||||
}
|
||||
} else {
|
||||
values3.push_back(-values2[j]);
|
||||
indices3.push_back(indices2[j]);
|
||||
}
|
||||
j++;
|
||||
}
|
||||
if (i<num1) {
|
||||
values3.insert(values3.end(), values1+i, values1+num1);
|
||||
indices3.insert(indices3.end(), indices1+i, indices1+num1);
|
||||
}
|
||||
if (j<num2) {
|
||||
for(index_t k=j; k<num2; k++) {
|
||||
values3.push_back(-values2[k]);
|
||||
indices3.push_back(indices2[k]);
|
||||
}
|
||||
}
|
||||
result->LoadRow(r, indices3, values3);
|
||||
}
|
||||
}
|
||||
/* Multiplication of two matrices A*B in matlab notation
|
||||
* If B is symmetric then it is much faster, because we can
|
||||
* multiply rows. Otherwise we have to compute the transpose
|
||||
* As an advise multiplication of two sparse matrices might
|
||||
* lead to a dense one, so please be carefull
|
||||
*/
|
||||
static inline void Multiply(const SparseMatrix &a,
|
||||
const SparseMatrix &b,
|
||||
SparseMatrix *result) {
|
||||
DEBUG_ASSERT(a.num_of_columns_ == b.num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_rows_ == result->num_of_rows_);
|
||||
DEBUG_ASSERT(b.num_of_columns_ == result->num_of_columns_);
|
||||
// If you try assigning the results to an already initialized matrix
|
||||
// you might get unexpected results. The following assertions
|
||||
// prevent you partially from that
|
||||
DEBUG_ASSERT(&a != result);
|
||||
DEBUG_ASSERT(&b != result);
|
||||
|
||||
if (b.issymmetric_ == true) {
|
||||
for(index_t r1=0; r1<a.num_of_rows_; r1++) {
|
||||
std::vector<index_t> indices3;
|
||||
std::vector<double> values3;
|
||||
index_t num1;
|
||||
double *values1;
|
||||
index_t *indices1;
|
||||
a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r1],
|
||||
num1, values1, indices1);
|
||||
for(index_t r2=0; r2<b.num_of_rows_; r2++) {
|
||||
index_t num2;
|
||||
double *values2;
|
||||
index_t *indices2;
|
||||
b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r2],
|
||||
num2, values2, indices2);
|
||||
index_t i=0;
|
||||
index_t j=0;
|
||||
double dot_product=0;
|
||||
while (likely(i<num1 && j<num2)) {
|
||||
while (indices1[i] < indices2[j]) {
|
||||
i++;
|
||||
if unlikely((i>=num1)) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (likely(i<num1) && indices1[i] == indices2[j]) {
|
||||
dot_product += values1[i] * values2[j];
|
||||
}
|
||||
j++;
|
||||
}
|
||||
if (dot_product!=0) {
|
||||
indices3.push_back(r2);
|
||||
values3.push_back(dot_product);
|
||||
}
|
||||
}
|
||||
result->LoadRow(r1, indices3, values3);
|
||||
indices3.clear();
|
||||
values3.clear();
|
||||
}
|
||||
} else {
|
||||
for(index_t r1=0; r1<a.num_of_rows_; r1++) {
|
||||
index_t num1;
|
||||
double *values1;
|
||||
index_t *indices1;
|
||||
a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r1],
|
||||
num1, values1, indices1);
|
||||
double dot_product=0;
|
||||
std::vector<index_t> indices3;
|
||||
std::vector<double> values3;
|
||||
for(index_t r2=0; r2<b.num_of_columns_; r2++) {
|
||||
for(index_t k=0; k< num1; k++) {
|
||||
dot_product += values1[k]*b.get(indices1[k], r2);
|
||||
}
|
||||
if (dot_product!=0){
|
||||
indices3.push_back(r2);
|
||||
values3.push_back(dot_product);
|
||||
}
|
||||
dot_product=0;
|
||||
}
|
||||
if (!indices3.empty()) {
|
||||
result->LoadRow(r1, indices3, values3);
|
||||
}
|
||||
indices3.clear();
|
||||
values3.clear();
|
||||
}
|
||||
}
|
||||
}
|
||||
/* The transpose flag should be set to true if
|
||||
* we want to use the transpose of mat, otherwise
|
||||
* set it to false.
|
||||
* */
|
||||
static inline void Multiply(const SparseMatrix &mat,
|
||||
const Vector &vec,
|
||||
Vector *result,
|
||||
bool transpose_flag) {
|
||||
Epetra_Vector temp_in(View, *(mat.map_), (double *)vec.ptr());
|
||||
Epetra_Vector temp_out(View, *(mat.map_), (double *)result->ptr());
|
||||
mat.matrix_->Multiply(transpose_flag, temp_in, temp_out);
|
||||
}
|
||||
/* Multiply the matrix with a scalar
|
||||
*/
|
||||
static inline void Multiply(const SparseMatrix &mat,
|
||||
const double scalar,
|
||||
SparseMatrix *result) {
|
||||
result->Copy(mat);
|
||||
result->Scale(scalar);
|
||||
}
|
||||
/* element wise multiplication of the matrices
|
||||
* A.*B in matlab notation
|
||||
*/
|
||||
static inline void DotMultiply(const SparseMatrix &a,
|
||||
const SparseMatrix &b,
|
||||
SparseMatrix *result) {
|
||||
DEBUG_ASSERT(a.num_of_columns_ == b.num_of_rows_);
|
||||
DEBUG_ASSERT(a.num_of_rows_ == result->num_of_rows_);
|
||||
DEBUG_ASSERT(b.num_of_columns_ == result->num_of_columns_);
|
||||
// If you try assigning the results to an already initialized matrix
|
||||
// you might get unexpected results. The following assertions
|
||||
// prevent you partially from that
|
||||
DEBUG_ASSERT(&a != result);
|
||||
DEBUG_ASSERT(&b != result);
|
||||
for(index_t r=0; r<a.num_of_rows_; r++) {
|
||||
std::vector<index_t> indices3;
|
||||
std::vector<double> values3;
|
||||
indices3.clear();
|
||||
values3.clear();
|
||||
index_t num1, num2;
|
||||
double *values1, *values2;
|
||||
index_t *indices1, *indices2;
|
||||
a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1);
|
||||
b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2);
|
||||
index_t i=0;
|
||||
index_t j=0;
|
||||
while (likely(i<num1 && j<num2)) {
|
||||
while (indices1[i] < indices2[j]) {
|
||||
i++;
|
||||
if unlikely((i>=num1)) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
if ( likely(i<num1) && indices1[i] == indices2[j]) {
|
||||
values3.push_back(values1[i] * values2[j]);
|
||||
indices3.push_back(indices1[i]);
|
||||
}
|
||||
j++;
|
||||
}
|
||||
result->LoadRow(r, indices3, values3);
|
||||
}
|
||||
}
|
||||
};
|
||||
#include "u/nvasil/sparse_matrix/sparse_matrix_impl.h"
|
||||
#endif
|
||||
@@ -0,0 +1,449 @@
|
||||
/*
|
||||
* =====================================================================================
|
||||
*
|
||||
* Filename: sparse_matrix_impl.h
|
||||
*
|
||||
* Description:
|
||||
*
|
||||
* Version: 1.0
|
||||
* Created: 12/02/2007 10:18:02 AM EST
|
||||
* Revision: none
|
||||
* Compiler: gcc
|
||||
*
|
||||
* Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org
|
||||
* Company: Georgia Tech Fastlab-ESP Lab
|
||||
*
|
||||
* =====================================================================================
|
||||
*
|
||||
*/
|
||||
|
||||
SparseMatrix::SparseMatrix() {
|
||||
map_ = NULL;
|
||||
issymmetric_=false;
|
||||
}
|
||||
|
||||
SparseMatrix::SparseMatrix(const index_t num_of_rows,
|
||||
const index_t num_of_columns,
|
||||
const index_t nnz_per_row) {
|
||||
Init(num_of_rows, num_of_columns, nnz_per_row);
|
||||
}
|
||||
void SparseMatrix::Init(const index_t num_of_rows,
|
||||
const index_t num_of_columns,
|
||||
const index_t nnz_per_row) {
|
||||
issymmetric_ = false;
|
||||
if likely(num_of_rows < num_of_columns) {
|
||||
FATAL("Num of rows %i should be greater than the num of columns %i\n",
|
||||
num_of_rows, num_of_columns);
|
||||
}
|
||||
num_of_rows_=num_of_rows;
|
||||
num_of_columns_=num_of_columns;
|
||||
dimension_ = num_of_rows_;
|
||||
if (num_of_rows_ < num_of_columns_) {
|
||||
FATAL("Num of rows %i is less than the number or columns %i",
|
||||
num_of_rows_, num_of_columns_);
|
||||
}
|
||||
map_ = new Epetra_Map(num_of_rows_, 0, comm_);
|
||||
matrix_ = Teuchos::rcp(
|
||||
new Epetra_CrsMatrix((Epetra_DataAccess)0, *map_, nnz_per_row));
|
||||
StartLoadingRows();
|
||||
}
|
||||
|
||||
void SparseMatrix::Init(index_t num_of_rows,
|
||||
index_t num_of_columns,
|
||||
index_t *nnz_per_row) {
|
||||
issymmetric_ = false;
|
||||
num_of_rows_=num_of_rows;
|
||||
num_of_columns_=num_of_columns;
|
||||
dimension_ = num_of_rows;
|
||||
if (num_of_rows_ < num_of_columns_) {
|
||||
FATAL("Num of rows %i is less than the number or columns %i",
|
||||
num_of_rows_, num_of_columns_);
|
||||
}
|
||||
map_ = new Epetra_Map(num_of_rows_, 0, comm_);
|
||||
matrix_ = Teuchos::rcp(
|
||||
new Epetra_CrsMatrix((Epetra_DataAccess)0 , *map_, nnz_per_row));
|
||||
StartLoadingRows();
|
||||
}
|
||||
|
||||
|
||||
void SparseMatrix::Init(const std::vector<index_t> &rows,
|
||||
const std::vector<index_t> &columns,
|
||||
const Vector &values,
|
||||
index_t nnz_per_row,
|
||||
index_t dimension) {
|
||||
issymmetric_ = false;
|
||||
if (nnz_per_row > 0 && dimension > 0) {
|
||||
Init(dimension, dimension, nnz_per_row);
|
||||
} else {
|
||||
num_of_rows_ = 0;
|
||||
num_of_columns_=0;
|
||||
std::map<index_t, index_t> frequencies;
|
||||
for(index_t i=0; i<(index_t)rows.size(); i++) {
|
||||
frequencies[rows[i]]++;
|
||||
frequencies[columns[i]]++;
|
||||
if (rows[i]>num_of_rows_) {
|
||||
num_of_rows_ = rows[i];
|
||||
}
|
||||
if (columns[i]>num_of_columns_) {
|
||||
num_of_columns_ = columns[i];
|
||||
}
|
||||
}
|
||||
num_of_columns_++;
|
||||
num_of_rows_++;
|
||||
if (num_of_rows_ < num_of_columns_) {
|
||||
FATAL("At this point we only support rows (%i) >= columns (%i)\n",
|
||||
num_of_rows_, num_of_columns_);
|
||||
}
|
||||
dimension_ = num_of_rows_;
|
||||
if ((index_t)frequencies.size()!=num_of_rows_) {
|
||||
NONFATAL("Some of the rows are zeros only!");
|
||||
}
|
||||
index_t *nnz= new index_t[dimension_];
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
nnz[i]=frequencies[i];
|
||||
}
|
||||
Init(num_of_rows_, num_of_columns_, nnz);
|
||||
//delete []nnz;
|
||||
}
|
||||
Load(rows, columns, values);
|
||||
}
|
||||
|
||||
void SparseMatrix::Init(const std::vector<index_t> &rows,
|
||||
const std::vector<index_t> &columns,
|
||||
const std::vector<double> &values,
|
||||
index_t nnz_per_row,
|
||||
index_t dimension) {
|
||||
issymmetric_ = false;
|
||||
Vector temp;
|
||||
temp.Alias((double *)&values[0], values.size());
|
||||
Init(rows, columns, temp, nnz_per_row, dimension);
|
||||
|
||||
}
|
||||
|
||||
void SparseMatrix::Init(std::string filename) {
|
||||
issymmetric_ = false;
|
||||
FILE *fp = fopen(filename.c_str(), "r");
|
||||
if (fp == NULL) {
|
||||
FATAL("Cannot open %s, error: %s",
|
||||
filename.c_str(),
|
||||
strerror(errno));
|
||||
}
|
||||
std::vector<index_t> rows;
|
||||
std::vector<index_t> cols;
|
||||
std::vector<double> vals;
|
||||
while (!feof(fp)) {
|
||||
index_t r, c;
|
||||
double v;
|
||||
fscanf(fp,"%i %i %lg\n", &r, &c, &v);
|
||||
rows.push_back(r);
|
||||
cols.push_back(c);
|
||||
vals.push_back(v);
|
||||
}
|
||||
fclose(fp);
|
||||
/*for(index_t i=0; i< (index_t)rows.size(); i++) {
|
||||
printf("%i %i %lg\n", rows[i], cols[i], vals[i]);
|
||||
}*/
|
||||
Init(rows, cols, vals, -1, -1);
|
||||
}
|
||||
|
||||
void SparseMatrix::Destruct() {
|
||||
if (map_!=NULL) {
|
||||
delete map_;
|
||||
map_=NULL;
|
||||
}
|
||||
}
|
||||
|
||||
void SparseMatrix::Copy(const SparseMatrix &other) {
|
||||
num_of_rows_ = other.num_of_rows_;
|
||||
num_of_columns_ = other.num_of_columns_;
|
||||
dimension_ = other.dimension_;
|
||||
map_ = new Epetra_Map(num_of_rows_, 0, comm_);
|
||||
issymmetric_ = other.issymmetric_;
|
||||
if (other.matrix_->Filled()==false) {
|
||||
matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(Epetra_DataAccess(0), *map_, 10));
|
||||
this->StartLoadingRows();
|
||||
for(index_t r=0; r<num_of_rows_; r++){
|
||||
index_t global_row = other.my_global_elements_[r];
|
||||
index_t num_of_entries;
|
||||
double *values;
|
||||
index_t *indices;
|
||||
other.matrix_->ExtractGlobalRowView(global_row, num_of_entries, values, indices);
|
||||
this->LoadRow(r, num_of_entries, indices, values);
|
||||
}
|
||||
} else {
|
||||
matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(*(other.matrix_.get())));
|
||||
}
|
||||
}
|
||||
void SparseMatrix::StartLoadingRows() {
|
||||
my_global_elements_ = map_->MyGlobalElements();
|
||||
}
|
||||
|
||||
void SparseMatrix::LoadRow(index_t row,
|
||||
std::vector<index_t> &columns,
|
||||
Vector &values) {
|
||||
DEBUG_ASSERT(values.length() == (index_t)columns.size());
|
||||
matrix_->InsertGlobalValues(my_global_elements_[row],
|
||||
values.length(),
|
||||
values.ptr(),
|
||||
&columns[0]);
|
||||
|
||||
}
|
||||
|
||||
void SparseMatrix::LoadRow(index_t row,
|
||||
index_t *columns,
|
||||
Vector &values) {
|
||||
matrix_->InsertGlobalValues(my_global_elements_[row],
|
||||
values.length(),
|
||||
values.ptr(),
|
||||
&columns[0]);
|
||||
|
||||
}
|
||||
|
||||
void SparseMatrix::LoadRow(index_t row,
|
||||
std::vector<index_t> &columns,
|
||||
std::vector<double> &values) {
|
||||
matrix_->InsertGlobalValues(my_global_elements_[row],
|
||||
values.size(),
|
||||
&values[0],
|
||||
&columns[0]);
|
||||
|
||||
}
|
||||
|
||||
void SparseMatrix::LoadRow(index_t row,
|
||||
index_t num,
|
||||
index_t *columns,
|
||||
double *values) {
|
||||
matrix_->InsertGlobalValues(my_global_elements_[row],
|
||||
num,
|
||||
values,
|
||||
columns);
|
||||
|
||||
}
|
||||
|
||||
void SparseMatrix::EndLoading() {
|
||||
matrix_->FillComplete();
|
||||
matrix_->OptimizeStorage();
|
||||
}
|
||||
|
||||
void SparseMatrix::Load(const std::vector<index_t> &rows,
|
||||
const std::vector<index_t> &columns,
|
||||
const Vector &values) {
|
||||
DEBUG_ASSERT(rows.size() ==columns.size());
|
||||
DEBUG_ASSERT((index_t)columns.size() == values.length());
|
||||
my_global_elements_ = map_->MyGlobalElements();
|
||||
index_t i=0;
|
||||
index_t cur_row = rows[i];
|
||||
index_t prev_row= rows[i];
|
||||
std::vector<index_t> indices;
|
||||
std::vector<double> row_values;
|
||||
while (true) {
|
||||
indices.clear();
|
||||
row_values.clear();
|
||||
while (likely((rows[cur_row]==rows[prev_row]) &&
|
||||
(i < (index_t)rows.size()))) {
|
||||
indices.push_back(columns[i]);
|
||||
row_values.push_back(values[i]);
|
||||
i++;
|
||||
prev_row=i-1;
|
||||
cur_row=i;
|
||||
}
|
||||
matrix_->InsertGlobalValues(my_global_elements_[rows[prev_row]],
|
||||
row_values.size(),
|
||||
&row_values[0],
|
||||
&indices[0]);
|
||||
prev_row=cur_row;
|
||||
if (i >= (index_t)rows.size()) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double SparseMatrix::get(index_t r, index_t c) const {
|
||||
DEBUG_BOUNDS(r, num_of_rows_);
|
||||
DEBUG_BOUNDS(c, num_of_columns_);
|
||||
index_t global_row = my_global_elements_[r];
|
||||
index_t num_of_entries;
|
||||
double *values;
|
||||
index_t *indices;
|
||||
matrix_->ExtractGlobalRowView(global_row, num_of_entries, values, indices);
|
||||
index_t *pos = std::find(indices, indices+num_of_entries, c);
|
||||
if (pos==indices+num_of_entries) {
|
||||
return 0;
|
||||
}
|
||||
return values[(ptrdiff_t)(pos-indices)];
|
||||
}
|
||||
|
||||
void SparseMatrix::set(index_t r, index_t c, double v) {
|
||||
DEBUG_BOUNDS(r, num_of_rows_);
|
||||
DEBUG_BOUNDS(c, num_of_columns_);
|
||||
if (get(r,c)!=0) {
|
||||
matrix_->InsertGlobalValues(my_global_elements_[r], 1, &v, &c);
|
||||
} else {
|
||||
matrix_->ReplaceGlobalValues(my_global_elements_[r], 1, &v, &c);
|
||||
}
|
||||
}
|
||||
void SparseMatrix::MakeSymmetric() {
|
||||
index_t num_of_entries;
|
||||
double *values;
|
||||
index_t *indices;
|
||||
for(index_t i=0; i<dimension_; i++) {
|
||||
index_t global_row = my_global_elements_[i];
|
||||
matrix_->ExtractGlobalRowView(global_row, num_of_entries, values, indices);
|
||||
for(index_t j=0; j<num_of_entries; j++) {
|
||||
if (unlikely(get(i, indices[j])!=values[j])) {
|
||||
set(i, indices[j], values[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
issymmetric_ = true;
|
||||
}
|
||||
|
||||
void SparseMatrix::Eig(index_t num_of_eigvalues,
|
||||
std::string eigtype,
|
||||
Matrix *eigvectors,
|
||||
std::vector<double> *real_eigvalues,
|
||||
std::vector<double> *imag_eigvalues) {
|
||||
if (unlikely(!matrix_->Filled())) {
|
||||
FATAL("You have to call EndLoading before running eigenvalues otherwise "
|
||||
"it will fail\n");
|
||||
}
|
||||
index_t block_size=2;
|
||||
Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp(new
|
||||
Epetra_MultiVector(*map_,
|
||||
block_size));
|
||||
// Fill it with random numbers
|
||||
ivec->Random();
|
||||
// Setup the eigenproblem, with the matrix A and the initial vectors ivec
|
||||
Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> >problem =
|
||||
Teuchos::rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(matrix_, ivec));
|
||||
|
||||
// The 2-D laplacian is symmetric. Specify this in the eigenproblem.
|
||||
if (issymmetric_ == true) {
|
||||
problem->setHermitian(true);
|
||||
} else {
|
||||
problem->setHermitian(false);
|
||||
}
|
||||
|
||||
// Specify the desired number of eigenvalues
|
||||
problem->setNEV(num_of_eigvalues);
|
||||
|
||||
// Signal that we are done setting up the eigenvalue problem
|
||||
bool ierr = problem->setProblem();
|
||||
|
||||
// Check the return from setProblem(). If this is true, there was an
|
||||
// error. This probably means we did not specify enough information for
|
||||
// the eigenproblem.
|
||||
if unlikely(ierr == false) {
|
||||
FATAL("Trilinos solver error, you probably didn't specify enough information "
|
||||
"for the eigenvalue problem\n");
|
||||
}
|
||||
|
||||
// Specify the verbosity level. Options include:
|
||||
// Anasazi::Errors
|
||||
// This option is always set
|
||||
// Anasazi::Warnings
|
||||
// Warnings (less severe than errors)
|
||||
// Anasazi::IterationDetails
|
||||
// Details at each iteration, such as the current eigenvalues
|
||||
// Anasazi::OrthoDetails
|
||||
// Details about orthogonality
|
||||
// Anasazi::TimingDetails
|
||||
// A summary of the timing info for the solve() routine
|
||||
// Anasazi::FinalSummary
|
||||
// A final summary
|
||||
// Anasazi::Debug
|
||||
// Debugging information
|
||||
int verbosity = Anasazi::Warnings +
|
||||
Anasazi::Errors +
|
||||
Anasazi::FinalSummary +
|
||||
Anasazi::TimingDetails;
|
||||
|
||||
// Choose which eigenvalues to compute
|
||||
// Choices are:
|
||||
// LM - target the largest magnitude [default]
|
||||
// SM - target the smallest magnitude
|
||||
// LR - target the largest real
|
||||
// SR - target the smallest real
|
||||
// LI - target the largest imaginary
|
||||
// SI - target the smallest imaginary
|
||||
|
||||
// Create the parameter list for the eigensolver
|
||||
Teuchos::ParameterList my_pl;
|
||||
my_pl.set( "Verbosity", verbosity);
|
||||
my_pl.set( "Which", eigtype);
|
||||
my_pl.set( "Block Size", block_size);
|
||||
my_pl.set( "Num Blocks", 20);
|
||||
my_pl.set( "Maximum Restarts", 100);
|
||||
my_pl.set( "Convergence Tolerance", 1.0e-8);
|
||||
|
||||
// Create the Block Krylov Schur solver
|
||||
// This takes as inputs the eigenvalue problem and the solver parameters
|
||||
Anasazi::BlockKrylovSchurSolMgr<double,MV,OP>
|
||||
my_block_krylov_schur(problem, my_pl);
|
||||
|
||||
// Solve the eigenvalue problem, and save the return code
|
||||
Anasazi::ReturnType solver_return = my_block_krylov_schur.solve();
|
||||
|
||||
// Check return code of the solver: Unconverged, Failed, or OK
|
||||
switch (solver_return) {
|
||||
// UNCONVERGED
|
||||
case Anasazi::Unconverged:
|
||||
NONFATAL("Anasazi::BlockKrylovSchur::solve() did not converge!\n");
|
||||
return ;
|
||||
// CONVERGED
|
||||
case Anasazi::Converged:
|
||||
NONFATAL("Anasazi::BlockKrylovSchur::solve() converged!\n");
|
||||
}
|
||||
// Get eigensolution struct
|
||||
Anasazi::Eigensolution<double, Epetra_MultiVector> sol = problem->getSolution();
|
||||
// Get the number of eigenpairs returned
|
||||
int num_of_eigvals_returned = sol.numVecs;
|
||||
if (num_of_eigvals_returned < num_of_eigvalues) {
|
||||
NONFATAL("The solver returned less eigenvalues (%i) "
|
||||
"than requested (%i)\n", num_of_eigvalues,
|
||||
num_of_eigvals_returned);
|
||||
}
|
||||
|
||||
// Get eigenvectors
|
||||
eigvectors->Init(dimension_, num_of_eigvals_returned);
|
||||
Teuchos::RCP<Epetra_MultiVector> evecs = sol.Evecs;
|
||||
evecs->ExtractCopy(eigvectors->GetColumnPtr(0), dimension_);
|
||||
|
||||
// Get eigenvalues
|
||||
std::vector<Anasazi::Value<double> > evals = sol.Evals;
|
||||
real_eigvalues->resize(num_of_eigvals_returned);
|
||||
for(index_t i=0; i<num_of_eigvals_returned; i++) {
|
||||
real_eigvalues->assign(i, evals[i].realpart);
|
||||
if (issymmetric_ == false) {
|
||||
imag_eigvalues->resize(num_of_eigvals_returned);
|
||||
imag_eigvalues->assign(i, evals[i].imagpart);
|
||||
}
|
||||
}
|
||||
|
||||
// Test residuals
|
||||
// Generate a (numev x numev) dense matrix for the eigenvalues...
|
||||
// This matrix is automatically initialized to zero
|
||||
Teuchos::SerialDenseMatrix<int, double> d(num_of_eigvalues,
|
||||
num_of_eigvalues);
|
||||
|
||||
// Add the eigenvalues on the diagonals (only the real part since problem is Hermitian)
|
||||
for (int i=0; i<num_of_eigvalues; i++) {
|
||||
d(i,i) = evals[i].realpart;
|
||||
}
|
||||
|
||||
// Generate a multivector for the product of the matrix and the eigenvectors
|
||||
Epetra_MultiVector res(*map_, num_of_eigvalues);
|
||||
|
||||
// R = A*evecs
|
||||
matrix_->Apply( *evecs, res);
|
||||
|
||||
// R -= evecs*D
|
||||
// = A*evecs - evecs*D
|
||||
MVT::MvTimesMatAddMv( -1.0, *evecs, d, 1.0, res);
|
||||
|
||||
// Compute the 2-norm of each vector in the MultiVector
|
||||
// and store them to a std::vector<double>
|
||||
std::vector<double> norm_res(num_of_eigvalues);
|
||||
MVT::MvNorm(res, &norm_res);
|
||||
}
|
||||
|
||||
@@ -0,0 +1,270 @@
|
||||
/*
|
||||
* =====================================================================================
|
||||
*
|
||||
* Filename: sparse_matrix_test.cc
|
||||
*
|
||||
* Description:
|
||||
*
|
||||
* Version: 1.0
|
||||
* Created: 12/08/2007 03:50:44 PM EST
|
||||
* Revision: none
|
||||
* Compiler: gcc
|
||||
*
|
||||
* Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org
|
||||
* Company: Georgia Tech Fastlab-ESP Lab
|
||||
*
|
||||
* =====================================================================================
|
||||
*/
|
||||
|
||||
#include <unistd.h>
|
||||
#include <errno.h>
|
||||
#include <sys/mman.h>
|
||||
#include <stdio.h>
|
||||
#include <limits>
|
||||
#include <vector>
|
||||
#include <map>
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "base/test.h"
|
||||
#include "sparse/sparse_matrix.h"
|
||||
|
||||
class SparseMatrixTest {
|
||||
public:
|
||||
SparseMatrixTest() {
|
||||
}
|
||||
~SparseMatrixTest() {
|
||||
}
|
||||
void Init() {
|
||||
for(index_t i=0; i<num_of_cols_; i++) {
|
||||
for(index_t j=0; j<num_of_rows_; j++) {
|
||||
mat_[i][j] = (i+j) * ((i+j) % 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
void Destruct() {
|
||||
delete smat_;
|
||||
}
|
||||
void TestInit1() {
|
||||
smat_ = new SparseMatrix(num_of_rows_,
|
||||
num_of_cols_,
|
||||
num_of_nnz_);
|
||||
smat_->StartLoadingRows();
|
||||
std::vector<index_t> ind;
|
||||
std::vector<double> val;
|
||||
for(index_t i=0; i<num_of_cols_; i++) {
|
||||
ind.clear();
|
||||
val.clear();
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
if (mat_[i][j] != 0) {
|
||||
ind.push_back(j);
|
||||
val.push_back(mat_[i][j]);
|
||||
}
|
||||
}
|
||||
smat_->LoadRow(i, ind, val);
|
||||
}
|
||||
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
TEST_DOUBLE_APPROX(smat_->get(i,j),
|
||||
mat_[i][j],
|
||||
std::numeric_limits<double>::epsilon());
|
||||
}
|
||||
}
|
||||
NONFATAL("TestInit1 sucess!!\n");
|
||||
}
|
||||
void TestInit2() {
|
||||
std::vector<index_t> rows;
|
||||
std::vector<index_t> cols;
|
||||
std::vector<double> vals;
|
||||
std::vector<index_t> nnz(num_of_rows_);
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
if (mat_[i][j] != 0) {
|
||||
rows.push_back(i);
|
||||
cols.push_back(j);
|
||||
vals.push_back(mat_[i][j]);
|
||||
nnz[i]++;
|
||||
}
|
||||
}
|
||||
}
|
||||
/*for(index_t i=0; i<(index_t)rows.size(); i++) {
|
||||
printf("%i %i %lg\n",
|
||||
rows[i],
|
||||
cols[i],
|
||||
vals[i]);
|
||||
}*/
|
||||
smat_ = new SparseMatrix();
|
||||
smat_->Init(rows, cols, vals,
|
||||
*(std::max_element(nnz.begin(), nnz.end())), num_of_rows_);
|
||||
// printf("%s\n", smat_->Print().c_str());
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
TEST_DOUBLE_APPROX(smat_->get(i,j), mat_[i][j],
|
||||
std::numeric_limits<double>::epsilon());
|
||||
}
|
||||
}
|
||||
NONFATAL("TestInit2 success!!\n");
|
||||
}
|
||||
void TestInit3() {
|
||||
FILE *fp = fopen("temp.txt", "w");
|
||||
if (fp==NULL) {
|
||||
FATAL("Cannot open temp.txt error %s", strerror(errno));
|
||||
}
|
||||
for(index_t i=0; i<num_of_cols_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
if (mat_[i][j] != 0) {
|
||||
fprintf(fp, "%i %i %g\n", i, j, mat_[i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
fclose(fp);
|
||||
smat_ = new SparseMatrix();
|
||||
smat_->Init("temp.txt");
|
||||
unlink("temp.txt");
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
TEST_DOUBLE_APPROX(smat_->get(i,j),
|
||||
mat_[i][j],
|
||||
std::numeric_limits<double>::epsilon());
|
||||
}
|
||||
}
|
||||
NONFATAL("TestInit3 success!!");
|
||||
}
|
||||
void TestCopyConstructor() {
|
||||
smat_ = new SparseMatrix();
|
||||
NONFATAL("TestCopyConstructor success!!\n");
|
||||
}
|
||||
void TestMakeSymmetric() {
|
||||
TestInit1();
|
||||
smat_->set(2, 3, 1.44);
|
||||
smat_->set(3, 2, 0.74);
|
||||
smat_->set(7, 8, 4.33);
|
||||
smat_->set(8, 7, 0.22);
|
||||
smat_->MakeSymmetric();
|
||||
for(index_t i=0; i<num_of_rows_; i++) {
|
||||
for(index_t j=0; j<num_of_cols_; j++) {
|
||||
TEST_DOUBLE_APPROX(smat_->get(i,j),
|
||||
smat_->get(j, i),
|
||||
std::numeric_limits<double>::epsilon());
|
||||
}
|
||||
}
|
||||
NONFATAL("Test MakeSymmetric success!!\n");
|
||||
}
|
||||
void TestEig() {
|
||||
TestInit1();
|
||||
smat_->EndLoading();
|
||||
std::vector<double> eigvalues_real;
|
||||
std::vector<double> eigvalues_imag;
|
||||
Matrix eigvectors;
|
||||
smat_->Eig(1, "LM", &eigvectors, &eigvalues_real, &eigvalues_imag);
|
||||
eigvectors.PrintDebug();
|
||||
NONFATAL("Test Eigenvector success!!\n");
|
||||
}
|
||||
void TestLinSolve() {
|
||||
TestInit1();
|
||||
Vector b,x;
|
||||
b.Init(num_of_cols_);
|
||||
b.SetZero();
|
||||
x.Init(num_of_cols_);
|
||||
x.SetAll(1);
|
||||
smat_->MakeSymmetric();
|
||||
smat_->EndLoading();
|
||||
smat_->LinSolve(b, &x);
|
||||
x.PrintDebug();
|
||||
NONFATAL("Test Linear Solve success!!\n");
|
||||
}
|
||||
void TestBasicOperations() {
|
||||
SparseMatrix a("A.txt");
|
||||
SparseMatrix b("B.txt");
|
||||
SparseMatrix a_plus_b("AplusB.txt");
|
||||
SparseMatrix a_minus_b("AminusB.txt");
|
||||
SparseMatrix a_times_b("AtimesB.txt");
|
||||
SparseMatrix a_dot_times_b("AdottimesB.txt");
|
||||
SparseMatrix temp;
|
||||
temp.Init(21,21, 3);
|
||||
Sparsem::Add(a, b, &temp);
|
||||
for(index_t i=0; i<20; i++) {
|
||||
for(index_t j=0; j<20; j++) {
|
||||
TEST_DOUBLE_APPROX(a_plus_b.get(i,j), temp.get(i,j), 0.01);
|
||||
}
|
||||
}
|
||||
temp.Destruct();
|
||||
NONFATAL("Matrix addition sucess!!\n");
|
||||
temp.Init(21,21, 3);
|
||||
Sparsem::Subtract(a, b, &temp);
|
||||
for(index_t i=0; i<21; i++) {
|
||||
for(index_t j=0; j<21; j++) {
|
||||
TEST_DOUBLE_APPROX(a_minus_b.get(i,j), temp.get(i,j), 0.01);
|
||||
}
|
||||
}
|
||||
temp.Destruct();
|
||||
NONFATAL("Matrix subtraction success!!\n");
|
||||
temp.Init(21,21, 3);
|
||||
Sparsem::Multiply(a, b, &temp);
|
||||
for(index_t i=0; i<21; i++) {
|
||||
for(index_t j=0; j<21; j++) {
|
||||
TEST_DOUBLE_APPROX(a_times_b.get(i,j), temp.get(i,j), 0.01);
|
||||
}
|
||||
}
|
||||
temp.Destruct();
|
||||
NONFATAL("Matrix multiplication success!!\n");
|
||||
temp.Init(21, 21, 3);
|
||||
Sparsem::DotMultiply(a, b, &temp);
|
||||
for(index_t i=0; i<a_dot_times_b.get_num_of_rows(); i++) {
|
||||
for(index_t j=0; j<a_dot_times_b.get_num_of_columns(); j++) {
|
||||
TEST_DOUBLE_APPROX(a_dot_times_b.get(i,j), temp.get(i,j), 0.01);
|
||||
}
|
||||
}
|
||||
temp.Destruct();
|
||||
NONFATAL("Matrix dot multiplication success!!\n");
|
||||
temp.Init(21,21, 3);
|
||||
Sparsem::Multiply(a, 3.45, &temp);
|
||||
for(index_t i=0; i<21; i++) {
|
||||
for(index_t j=0; j<21; j++) {
|
||||
TEST_DOUBLE_APPROX(3.45 * a.get(i,j), temp.get(i,j),
|
||||
std::numeric_limits<double>::epsilon());
|
||||
}
|
||||
}
|
||||
temp.Destruct();
|
||||
NONFATAL("Matrix scalar multiplicationn success!!\n");
|
||||
|
||||
}
|
||||
void TestAll() {
|
||||
Init();
|
||||
TestInit1();
|
||||
Destruct();
|
||||
Init();
|
||||
TestInit2();
|
||||
Destruct();
|
||||
Init();
|
||||
TestInit3();
|
||||
Destruct();
|
||||
Init();
|
||||
TestCopyConstructor();
|
||||
Destruct();
|
||||
Init();
|
||||
TestMakeSymmetric();
|
||||
Destruct();
|
||||
Init();
|
||||
TestEig();
|
||||
Destruct();
|
||||
Init();
|
||||
TestLinSolve();
|
||||
Destruct();
|
||||
Init();
|
||||
TestBasicOperations();
|
||||
}
|
||||
|
||||
private:
|
||||
SparseMatrix *smat_;
|
||||
static const index_t num_of_cols_ = 80;
|
||||
static const index_t num_of_rows_ = 80;
|
||||
static const index_t num_of_nnz_ = 4;
|
||||
double mat_[num_of_rows_][num_of_cols_];
|
||||
std::vector<index_t> indices_;
|
||||
std::vector<index_t> rows_;
|
||||
};
|
||||
|
||||
int main() {
|
||||
SparseMatrixTest test;
|
||||
test.TestAll();
|
||||
}
|
||||
Reference in New Issue
Block a user