Homec4science
Diffusion Eigenvalue Problems (master) Browse

eigpro/Matrix.cpp
a0a6218630c7master

Last Change
Raw File

Matrix.cpp
View Options

#include "Matrix.h"
#include <string>
#include <fstream>
#include <math.h>
#include <cstdlib>
#include <cmath>
#include <tgmath.h>
#include <iostream>
using namespace std;
Matrix::Matrix (string filename){
int x, y;
ifstream in(filename.c_str());
if (!in) {
cout << "Cannot open file.\n";
return;
}
in >> size;
matrix = new double*[size];
for (int i = 0; i < size; ++i)
matrix[i] = new double[size];
for (x = 0; x < size; x++) {
for (y = 0; y < size; y++) {
in >> matrix[x][y];
}
}
in >> shift_amount;
in.close();
}
Matrix::~Matrix(){
int i, j;
for(i = 0 ; i < size ; i++) {
delete[] matrix[i]; // Delete the rows of matrix
}
// Finally delete the matrix itself
delete[] matrix;
}
void Matrix::printMatrix(double** matrix, int size){
int x, y;
cout << endl << "Matrix Content:" << endl;
for(x = 0; x < size; x++) // loop 3 times for three lines
{
for(y = 0; y < size; y++) // loop for the three elements on the line
{
cout<<matrix[x][y]<<" "; // display the current element out of the array
}
cout<<endl; // when the inner loop is done, go to a new line
}
return;
}
double Matrix::powerMethod(double** m, int size, double tol, double shift_amount) {
int i, j;
double c[size], x[size], d, tmp;
double **matrix;
d = 0;
if(shift_amount != 0){
matrix = shift(m, size, shift_amount);
}
else {
matrix = m;
}
// We start with an initial approximation of the dominant eigenvector of A, which must be non-zero
for (i = 0; i < size; i++) {
x[i] = rand() % 3 + 1;
}
do
{
// c = A * x^k
for(i = 0; i < size; i++)
{
c[i] = 0;
for(j = 0; j < size; j++) {
c[i] += matrix[i][j] * x[j];
}
}
// x = c
for(i = 0; i < size; i++) {
x[i] = c[i];
}
// tmp = d_old
tmp = d;
// d = norm(x)
d = linf_norm(x, size);
// x = x/norm(x)
for(i = 0; i < size; i++) {
x[i] /= d;
}
}while(fabs(d-tmp) > tol);
cout << endl << "POWER METHOD:" << endl;
if(shift_amount == 0) {
cout << "Eigenvalue is: " << d << endl;
}
else{
cout << "Eigenvalue is: " << shift_amount + d << endl;
}
cout << "Eigenvector is: ";
for(i = 0; i < size; i++) {
cout << endl << x[i];
}
cout << endl;
if(shift_amount == 0) {
return d;
}
else {
return shift_amount + d;
}
}
double Matrix::inversePowerMethod(double** matrix, int size, double tol, double shift_amount){
int P[size+1], i;
double x[size], c[size], d, tmp;
d = 0;
if(!LUPDecompose(matrix, size, tol, P)){
cout << "Failure, Matrix is degenerate!" << endl;
exit(-1);
};
// We start with an initial approximation of the dominant eigenvector of A, which must be non-zero
for (i = 0; i < size; i++) {
x[i] = rand() % 3 + 1;
}
do
{
// c = A^(-1) * x^k
// A * c = x^k
LUPSolve(matrix, P, x, size, c);
// x = c
for(i = 0; i < size; i++) {
x[i] = c[i];
}
// tmp = d_old
tmp = d;
// d = norm(x)
d = linf_norm(x, size);
// x = x/norm(x)
for(i = 0; i < size; i++) {
x[i] /= d;
}
}while(fabs(d-tmp) > tol);
cout << endl << "INVERSE POWER METHOD:" << endl;
if(shift_amount == 0) {
cout << "Eigenvalue is: " << 1 / d << endl;
}
else{
cout << "Eigenvalue is: " << shift_amount + 1 / d << endl;
}
cout<<"Eigenvector is: ";
for(i = 0; i < size; i++) {
cout << endl << x[i];
}
cout << endl;
if(shift_amount == 0) {
return 1/d;
}
else {
return shift_amount + 1/d;
}
}
double Matrix::l2_norm(double *u, int size) {
double accum = 0.;
int i;
for (i = 0; i < size; i++) {
accum += u[i] * u[i];
}
return sqrt(accum);
}
double Matrix::linf_norm(double *u, int size) {
double d = 0.;
int i;
for(i = 0; i < size; i++)
{
if(fabs(u[i])>fabs(d))
d=u[i];
}
return d;
}
double** Matrix::shift(double** A, int size, double shift){
int i, j;
for(i = 0; i < size; i++)
{
for(j = 0; j < size; j++){
if(j == i){
A[i][j] -= shift;
}
}
}
return A;
}
int Matrix::LUPDecompose(double **A, int N, double Tol, int *P) {
int i, j, k, imax;
double maxA, *ptr, absA;
for (i = 0; i <= N; i++)
P[i] = i; //Unit permutation matrix, P[N] initialized with N
for (i = 0; i < N; i++) {
maxA = 0;
imax = i;
for (k = i; k < N; k++)
if ((absA = fabs(A[k][i])) > maxA) {
maxA = absA;
imax = k;
}
if (maxA < Tol) return 0; //failure, matrix is degenerate
if (imax != i) {
//pivoting P
j = P[i];
P[i] = P[imax];
P[imax] = j;
//pivoting rows of A
ptr = A[i];
A[i] = A[imax];
A[imax] = ptr;
//counting pivots starting from N (for determinant)
P[N]++;
}
for (j = i + 1; j < N; j++) {
A[j][i] /= A[i][i];
for (k = i + 1; k < N; k++)
A[j][k] -= A[j][i] * A[i][k];
}
}
return 1; //decomposition done
}
void Matrix::LUPSolve(double **A, int *P, double *b, int N, double *x) {
int i;
for (i = 0; i < N; i++) {
x[i] = b[P[i]];
for (int k = 0; k < i; k++)
x[i] -= A[i][k] * x[k];
}
for (i = N - 1; i >= 0; i--) {
for (int k = i + 1; k < N; k++)
x[i] -= A[i][k] * x[k];
x[i] = x[i] / A[i][i];
}
return;
}
int Matrix::get_size() {
return size;
}
double** Matrix::get_matrix() {
return matrix;
}
double Matrix::get_shift_amount() {
return shift_amount;
}
int Matrix::test(int i) {
double tol = 0.00000000001, eigenvalue;
// Compute eigenvalue using Power Method
eigenvalue = powerMethod(get_matrix(), get_size(), tol, get_shift_amount());
switch (i) {
case 1:
if (round(eigenvalue)==4){
cout << "Test 1: Success!" << endl;
}
else{
cout << "Test 1: Failure!" << endl;
}
break;
case 2:
if (round(eigenvalue)==2){
cout << "Test 2: Success!" << endl;
}
else{
cout << "Test 2: Failure!" << endl;
}
break;
case 3:
if (round(eigenvalue)==2){
cout << "Test 3: Success!" << endl;
}
else{
cout << "Test 3: Failure!" << endl;
}
break;
case 4:
if (round(eigenvalue)==2){
cout << "Test 4: Success!" << endl;
}
else{
cout << "Test 4: Failure!" << endl;
}
break;
case 5:
if (round(eigenvalue)==2){
cout << "Test 5: Success!" << endl;
}
else{
cout << "Test 5: Failure!" << endl;
}
break;
case 6:
if (round(eigenvalue)==3.8){
cout << "Test 6: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is 3.8" << endl;
}
break;
case 7:
if (round(eigenvalue)==3.8){
cout << "Test 7: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is 3.8" << endl;
}
break;
default:
break;
}
// Compute eigenvalue using Inverse Power Method
eigenvalue = inversePowerMethod(get_matrix(), get_size(), tol, get_shift_amount());
switch (i) {
case 1:
if (round(eigenvalue)==1){
cout << "Test 1: Success!" << endl;
}
else{
cout << "Test 1: Failure!" << endl;
}
break;
case 2:
if (round(eigenvalue)==1){
cout << "Test 2: Success!" << endl;
}
else{
cout << "Test 2: Failure!" << endl;
}
break;
case 3:
if (round(eigenvalue)==-1.5){
cout << "Test 3: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is -1.5" << endl;
}
break;
case 4:
if (round(eigenvalue)==1){
cout << "Test 4: Success!" << endl;
}
else{
cout << "Test 4: Failure!" << endl;
}
break;
case 5:
if (round(eigenvalue)==1){
cout << "Test 5: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is 1" << endl;
}
break;
case 6:
if (round(eigenvalue)==0.2){
cout << "Test 6: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is 0.2!" << endl;
}
break;
case 7:
if (round(eigenvalue)==0.75){
cout << "Test 7: Success!" << endl;
}
else{
cout << "The Eigenvalue calculated by WolframAlpha is 0.75" << endl;
}
break;
default:
break;
}
}
c4science · Help