Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88113963
matrixfun.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, Oct 16, 21:05
Size
8 KB
Mime Type
text/x-c
Expires
Fri, Oct 18, 21:05 (2 d)
Engine
blob
Format
Raw Data
Handle
21710882
Attached To
rLAMMPS lammps
matrixfun.cpp
View Options
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: matrixfun.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "matrixfun.h"
#include <math.h>
#include "fastmatrixops.h"
#include <cstdlib>
using namespace std;
//
// Create a new matrix
//
VirtualMatrix* NewMatrix(int type){
switch( MatrixType(type) )
{
case MATRIX : return new Matrix;
case COLMATRIX : return new ColMatrix;
case ROWMATRIX : return new RowMatrix;
case MAT3X3 : return new Mat3x3;
case MAT4X4 : return new Mat4x4;
case VECT3 : return new Vect3;
case VECT4 : return new Vect4;
default : return 0; // error!
}
}
//
// Transpose
//
Matrix T(const VirtualMatrix& A){
int numrows = A.GetNumRows();
int numcols = A.GetNumCols();
Matrix C(numcols,numrows);
for(int i=0;i<numcols;i++)
for(int j=0;j<numrows;j++)
C.BasicSet(i,j,A.BasicGet(j,i));
return C;
}
Mat3x3 T(const Mat3x3& A){
Mat3x3 C;
C.elements[0][0] = A.elements[0][0];
C.elements[1][1] = A.elements[1][1];
C.elements[2][2] = A.elements[2][2];
C.elements[0][1] = A.elements[1][0];
C.elements[0][2] = A.elements[2][0];
C.elements[1][2] = A.elements[2][1];
C.elements[1][0] = A.elements[0][1];
C.elements[2][0] = A.elements[0][2];
C.elements[2][1] = A.elements[1][2];
return C;
}
Mat6x6 T(const Mat6x6& A){
Mat6x6 C;
int i,j;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
C.elements[i][j] = A.elements[j][i];
return C;
}
Matrix T(const Vect3& A){
Matrix C(1,3);
C.BasicSet(0,0,A.elements[0]);
C.BasicSet(0,1,A.elements[1]);
C.BasicSet(0,2,A.elements[2]);
return C;
}
Matrix T(const Vect6& A){
Matrix C(1,6);
C.BasicSet(0,0,A.elements[0]);
C.BasicSet(0,1,A.elements[1]);
C.BasicSet(0,2,A.elements[2]);
C.BasicSet(0,3,A.elements[3]);
C.BasicSet(0,4,A.elements[4]);
C.BasicSet(0,5,A.elements[5]);
return C;
}
RowMatrix T(const VirtualColMatrix &A){
int numele = A.GetNumRows();
RowMatrix C(numele);
for(int i=0;i<numele;i++)
C.BasicSet(i,A.BasicGet(i));
return C;
}
ColMatrix T(const VirtualRowMatrix &A){
int numele = A.GetNumCols();
ColMatrix C(numele);
for(int i=0;i<numele;i++)
C.BasicSet(i,A.BasicGet(i));
return C;
}
//
// Symmetric Inverse
//
Matrix SymInverse(Matrix &A){
int r = A.GetNumRows();
Matrix C(r,r);
Matrix LD(r,r);
Matrix I(r,r);
I.Zeros();
for(int i=0;i<r;i++) I.BasicSet(i,i,1.0);
FastLDLT(A,LD);
FastLDLTSubs(LD,I,C);
return C;
}
Mat6x6 SymInverse(Mat6x6 &A){
Mat6x6 C;
Mat6x6 LD;
Mat6x6 I;
I.Zeros();
for(int i=0;i<6;i++) I.BasicSet(i,i,1.0);
FastLDLT(A,LD);
FastLDLTSubs(LD,I,C);
return C;
}
//
// Inverse
//
Matrix Inverse(Matrix& A){
int r = A.GetNumRows();
int indx[10000];
Matrix LU(r,r);
Matrix I(r,r);
Matrix C(r,r);
I.Zeros();
for(int i=0;i<r;i++) I.BasicSet(i,i,1.0);
FastLU(A,LU,indx);
FastLUSubs(LU,I,C,indx);
return C;
}
Mat3x3 Inverse(Mat3x3& A){
int indx[10000];
Mat3x3 LU;
Matrix I(3,3);
Matrix C(3,3);
I.Zeros();
for(int i=0;i<3;i++) I.BasicSet(i,i,1.0);
FastLU(A,LU,indx);
FastLUSubs(LU,I,C,indx);
return C;
}
Mat4x4 Inverse(Mat4x4& A){
int indx[10000];
Mat4x4 LU;
Matrix I(4,4);
Matrix C(4,4);
I.Zeros();
for(int i=0;i<4;i++) I.BasicSet(i,i,1.0);
FastLU(A,LU,indx);
FastLUSubs(LU,I,C,indx);
return C;
}
Mat6x6 Inverse(Mat6x6& A){
int indx[10000];
Mat6x6 LU;
Matrix I(6,6);
Matrix C(6,6);
I.Zeros();
for(int i=0;i<6;i++) I.BasicSet(i,i,1.0);
FastLU(A,LU,indx);
FastLUSubs(LU,I,C,indx);
return C;
}
//
// overloaded addition
//
Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B){ // addition
int Arows,Acols,Brows,Bcols;
Arows = A.GetNumRows();
Acols = A.GetNumCols();
Brows = B.GetNumRows();
Bcols = B.GetNumCols();
if( !((Arows == Brows) && (Acols == Bcols)) ){
cerr << "Dimesion mismatch in matrix addition" << endl;
exit(1);
}
Matrix C(Arows,Acols);
for(int i=0;i<Arows;i++)
for(int j=0;j<Acols;j++)
C.BasicSet(i,j,A.BasicGet(i,j) + B.BasicGet(i,j));
return C;
}
//
// overloaded subtraction
//
Matrix operator- (const VirtualMatrix &A, const VirtualMatrix &B){ // subtraction
int Arows,Acols,Brows,Bcols;
Arows = A.GetNumRows();
Acols = A.GetNumCols();
Brows = B.GetNumRows();
Bcols = B.GetNumCols();
if( !((Arows == Brows) && (Acols == Bcols)) ){
cerr << "Dimesion mismatch in matrix addition" << endl;
exit(1);
}
Matrix C(Arows,Acols);
for(int i=0;i<Arows;i++)
for(int j=0;j<Acols;j++)
C.BasicSet(i,j,A.BasicGet(i,j) - B.BasicGet(i,j));
return C;
}
//
// overloaded matrix multiplication
//
Matrix operator* (const VirtualMatrix &A, const VirtualMatrix &B){ // multiplication
int Arows,Acols,Brows,Bcols;
Arows = A.GetNumRows();
Acols = A.GetNumCols();
Brows = B.GetNumRows();
Bcols = B.GetNumCols();
if(Acols != Brows){
cerr << "Dimesion mismatch in matrix multiplication" << endl;
exit(1);
}
Matrix C(Arows,Bcols);
C.Zeros();
for(int i=0;i<Arows;i++)
for(int j=0;j<Bcols;j++)
for(int k=0;k<Brows;k++)
C.BasicIncrement(i,j, A.BasicGet(i,k) * B.BasicGet(k,j) );
return C;
}
//
// overloaded scalar multiplication
//
Matrix operator* (const VirtualMatrix &A, double b){ // multiplication
Matrix C = A;
C *= b;
return C;
}
Matrix operator* (double b, const VirtualMatrix &A){ // multiplication
Matrix C = A;
C *= b;
return C;
}
//
// overloaded negative
//
Matrix operator- (const VirtualMatrix& A){ // negative
int r = A.GetNumRows();
int c = A.GetNumCols();
Matrix C(r,c);
for(int i=0;i<r;i++)
for(int j=0;j<c;j++)
C.BasicSet(i,j,-A.BasicGet(i,j));
return C;
}
//
// Cross product (friend of Vect3)
//
Vect3 Cross(Vect3& a, Vect3& b){
return CrossMat(a)*b;
}
//
// Cross Matrix (friend of Vect3 & Mat3x3)
//
Mat3x3 CrossMat(Vect3& a){
Mat3x3 C;
C.Zeros();
C.elements[0][1] = -a.elements[2];
C.elements[1][0] = a.elements[2];
C.elements[0][2] = a.elements[1];
C.elements[2][0] = -a.elements[1];
C.elements[1][2] = -a.elements[0];
C.elements[2][1] = a.elements[0];
return C;
}
//
// Stack
//
Matrix Stack(VirtualMatrix& A, VirtualMatrix& B){
int m,na,nb;
m = A.GetNumCols();
if( m != B.GetNumCols()){
cerr << "Error: cannot stack matrices of differing column dimension" << endl;
exit(0);
}
na = A.GetNumRows();
nb = B.GetNumRows();
Matrix C(na+nb,m);
for(int i=0;i<na;i++){
for(int j=0;j<m;j++){
C.BasicSet(i,j,A.BasicGet(i,j));
}
}
for(int i=0;i<nb;i++){
for(int j=0;j<m;j++){
C.BasicSet(i+na,j,B.BasicGet(i,j));
}
}
return C;
}
//Hstack
Matrix HStack(VirtualMatrix& A, VirtualMatrix& B){
int m,na,nb;
m = A.GetNumRows();
if( m != B.GetNumRows()){
cerr << "Error: cannot stack matrices of differing row dimension" << endl;
exit(0);
}
na = A.GetNumCols();
nb = B.GetNumCols();
Matrix C(m,na+nb);
for(int i=0;i<m;i++){
for(int j=0;j<na;j++){
C.BasicSet(i,j,A.BasicGet(i,j));
}
}
for(int i=0;i<m;i++){
for(int j=0;j<nb;j++){
C.BasicSet(i,j+na,B.BasicGet(i,j));
}
}
return C;
}
//
//
void Set6DAngularVector(Vect6& v6, Vect3& v3){
v6.elements[0] = v3.elements[0];
v6.elements[1] = v3.elements[1];
v6.elements[2] = v3.elements[2];
}
void Set6DLinearVector(Vect6& v6, Vect3& v3){
v6.elements[3] = v3.elements[0];
v6.elements[4] = v3.elements[1];
v6.elements[5] = v3.elements[2];
}
Event Timeline
Log In to Comment