Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76377019
Matrix.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Aug 7, 15:30
Size
5 KB
Mime Type
text/x-c
Expires
Fri, Aug 9, 15:30 (2 d)
Engine
blob
Format
Raw Data
Handle
19584251
Attached To
rLAMMPS lammps
Matrix.cpp
View Options
#include "DenseMatrix.h"
#include "Solver.h"
//-----------------------------------------------------------------------------
//* performs a matrix-matrix multiply with optional transposes BLAS version
// C = b*C + a*A*B
//-----------------------------------------------------------------------------
void MultAB(const MATRIX &A, const MATRIX &B, DENS_MAT &C,
const bool At, const bool Bt, double a, double b)
{
static char t[2] = {'N','T'};
char *ta=t+At, *tb=t+Bt;
int sA[2] = {A.nRows(), A.nCols()}; // sizes of A
int sB[2] = {B.nRows(), B.nCols()}; // sizes of B
GCK(A, B, sA[!At]!=sB[Bt], "MultAB<double>: matrix-matrix multiply");
if (!C.is_size(sA[At],sB[!Bt]))
{
C.resize(sA[At],sB[!Bt]);
if (b != 0.0) C.zero();
}
// get pointers to the matrix sizes needed by BLAS
int *M = sA+At; // # of rows in op[A] (op[A] = A' if At='T' else A)
int *N = sB+!Bt; // # of cols in op[B]
int *K = sA+!At; // # of cols in op[A] or # of rows in op[B]
double *pa=A.get_ptr(), *pb=B.get_ptr(), *pc=C.get_ptr();
#ifdef COL_STORAGE
dgemm_(ta, tb, M, N, K, &a, pa, sA, pb, sB, &b, pc, M);
#else
dgemm_(tb, ta, N, M, K, &a, pb, sB+1, pa, sA+1, &b, pc, N);
#endif
}
//-----------------------------------------------------------------------------
//* returns the inverse of a double precision matrix
//-----------------------------------------------------------------------------
DenseMatrix<double> inv(const MATRIX& A)
{
SQCK(A, "DenseMatrix::inv(), matrix not square"); // check matrix is square
DENS_MAT invA(A); // Make copy of A to invert
// setup for call to BLAS
int m, info, lwork=-1;
m = invA.nRows();
int *ipiv = new int[m<<1]; // need 2m storage
int *iwork=ipiv+m;
dgetrf_(&m,&m,invA.get_ptr(),&m,ipiv,&info); // compute LU factorization
GCK(A,A,info<0,"DenseMatrix::inv() dgetrf error: Argument had bad value.");
GCK(A,A,info>0,"DenseMatrix::inv() dgetrf error: Matrix not invertable.");
// LU factorization succeeded
// Compute 1-norm of original matrix for use with dgecon
char norm = '1'; // Use 1-norm
double rcond, anorm, *workc = new double[4*m];
anorm = dlange_(&norm,&m,&m,A.get_ptr(),&m,workc);
// Condition number estimation (warn if bad)
dgecon_(&norm,&m,invA.get_ptr(),&m,&anorm,&rcond,workc,iwork,&info);
GCK(A,A,info<0, "DenseMatrix::inv(): dgecon error: Argument had bad value.");
GCK(A,A,rcond<1e-14,"DenseMatrix::inv(): Matrix nearly singular, RCOND<e-14");
// Now determine optimal work size for computation of matrix inverse
double work_dummy[2] = {0.0,0.0};
dgetri_(&m, invA.get_ptr(), &m, ipiv, work_dummy, &lwork, &info);
GCK(A,A,info<0,"DenseMatrix::inv() dgetri error: Argument had bad value.");
GCK(A,A,info>0,"DenseMatrix::inv() dgetri error: Matrix not invertable.");
// Work size query succeded
lwork = (int)work_dummy[0];
double *work = new double[lwork]; // Allocate vector of appropriate size
// Compute and store matrix inverse
dgetri_(&m,invA.get_ptr(),&m,ipiv,work,&lwork,&info);
GCK(A,A,info<0,"DenseMatrix::inv() dgetri error: Argument had bad value.");
GCK(A,A,info>0,"DenseMatrix::inv() dgetri error: Matrix not invertable.");
// Clean-up
delete [] ipiv;
delete [] workc;
delete [] work;
return invA;
}
//-----------------------------------------------------------------------------
//* computes the determinant of a square matrix by LU decomposition (if n>3)
//-----------------------------------------------------------------------------
double det(const MATRIX& A)
{
static const double sign[2] = {1.0, -1.0};
SQCK(A, "Matrix::det(), matrix not square"); // check matrix is square
int m = A.nRows();
switch (m) // explicit determinant for small matrix sizes
{
case 1: return A(0,0);
case 2: return A(0,0)*A(1,1)-A(0,1)*A(1,0);
case 3:
return A(0,0)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))
+ A(0,1)*(A(1,2)*A(2,0)-A(1,0)*A(2,2))
+ A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
default: break;
}
// First compute LU factorization
int info, *ipiv = new int[m];
double det = 1.0;
DENS_MAT PLUA(A);
dgetrf_(&m,&m,PLUA.get_ptr(),&m,ipiv,&info);
GCK(A,A,info>0,"Matrix::det() dgetrf error: Bad argument value.");
if (!info) // matrix is non-singular
{
// Compute det(A) = det(P)*det(L)*det(U) = +/-1 * det(U)
int i, OddNumPivots;
det = PLUA(0,0);
OddNumPivots = ipiv[0]!=(1);
for(i=1; i<m; i++)
{
det *= PLUA(i,i);
OddNumPivots += (ipiv[i]!=(i+1)); // # pivots even/odd
}
det *= sign[OddNumPivots&1];
}
delete [] ipiv; // Clean-up
return det;
}
//-----------------------------------------------------------------------------
//* Returns the maximum eigenvalue of a matrix.
//-----------------------------------------------------------------------------
double max_eigenvalue(const Matrix<double>& A)
{
GCK(A,A,!is_size(3,3), "max_eigenvalue only implimented for 3x3");
const double c0 = det(A);
const double c1 = A(1,0)*A(0,1) + A(2,0)*A(0,2) + A(1,2)*A(2,1)
- A(0,0)*A(1,1) - A(0,0)*A(2,2) - A(1,1)*A(2,2);
const double c2 = trace(A);
double c[4] = {c0, c1, c2, -1.0}, x[3];
int num_roots = ATC::solve_cubic(c, x);
double max_root = 0.0;
for (int i=0; i<num_roots; i++) max_root = std::max(x[i], max_root);
return max_root;
}
Event Timeline
Log In to Comment