You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1000 lines
25 KiB
1000 lines
25 KiB
#include <cstdio>
|
|
#include <cstdlib>
|
|
#include <string>
|
|
#include <iostream>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include "Constants.h"
|
|
#include "sparmat.h"
|
|
#include "densemat.h"
|
|
#include "hb_io.h"
|
|
|
|
using namespace std;
|
|
|
|
sparMat::sparMat()
|
|
{
|
|
type = SYMMETRIC;
|
|
dim = 0;
|
|
colDim = 0;
|
|
NZ = 0;
|
|
NUMSEG = 0;
|
|
SEGSIZE = 0;
|
|
rowA = NULL;
|
|
colA = NULL;
|
|
entry = NULL;
|
|
}
|
|
|
|
sparMat::~sparMat()
|
|
{
|
|
freeMemory();
|
|
}
|
|
|
|
void sparMat::freeMemory()
|
|
{
|
|
delete [] rowA; rowA = NULL;
|
|
delete [] colA; colA = NULL;
|
|
if (entry != NULL) {
|
|
delete [] entry[0];
|
|
// for (int i = 1; i < NUMSEG; i++)
|
|
// entry[i] = new fp_t;
|
|
delete [] entry; entry = NULL;
|
|
}
|
|
}
|
|
|
|
void sparMat::init(int n, int nonZero)
|
|
{
|
|
type = SYMMETRIC;
|
|
dim = n;
|
|
colDim = n;
|
|
NZ = nonZero;
|
|
NUMSEG = 1;
|
|
SEGSIZE = NZ/NUMSEG;
|
|
rowA = new int[dim + 1];
|
|
colA = new int[NZ];
|
|
entry = new fp_t*[NUMSEG];
|
|
entry[0] = new fp_t[NZ];
|
|
// for (int i = 1; i < NUMSEG; i++)
|
|
// entry[i] = new fp_t;
|
|
}
|
|
|
|
void sparMat::sparInit(int n, int nonZero)
|
|
{
|
|
type = SYMMETRIC;
|
|
dim = n;
|
|
colDim = n;
|
|
NZ = nonZero;
|
|
NUMSEG = 1;
|
|
SEGSIZE = NZ/NUMSEG;
|
|
rowA = new int[dim + 1];
|
|
colA = new int[NZ];
|
|
}
|
|
|
|
void sparMat::sparInit(int n, int m, int nonZero)
|
|
{
|
|
type = NON_SYMMETRIC;
|
|
dim = n;
|
|
colDim = m;
|
|
NZ = nonZero;
|
|
NUMSEG = 1;
|
|
SEGSIZE = NZ/NUMSEG;
|
|
entry = new fp_t*[NUMSEG];
|
|
entry[0] = new fp_t[NZ];
|
|
// for (int i = 1; i < NUMSEG; i++)
|
|
// entry[i] = new fp_t;
|
|
}
|
|
|
|
// debug
|
|
void sparMat::sparInitNonSymm(int n, int nonZero)
|
|
{
|
|
type = NON_SYMMETRIC;
|
|
dim = n;
|
|
colDim = n;
|
|
NZ = nonZero;
|
|
NUMSEG = 1;
|
|
SEGSIZE = NZ/NUMSEG;
|
|
rowA = new int[dim + 1];
|
|
colA = new int[NZ];
|
|
entry = NULL;
|
|
}
|
|
|
|
void sparMat::sparINIT(int n, int m, int nonZero)
|
|
{
|
|
type = NON_SYMMETRIC;
|
|
dim = n;
|
|
colDim = m;
|
|
NZ = nonZero;
|
|
rowA = new int[dim + 1];
|
|
colA = new int[NZ];
|
|
NUMSEG = 1;
|
|
SEGSIZE = NZ;
|
|
entry = new fp_t*[NUMSEG];
|
|
entry[0] = new fp_t[NZ]; // orignal
|
|
// Modified by Stelios
|
|
for (int i = 0; i < NZ; i++)
|
|
entry[0][i] = 0.0;
|
|
}
|
|
|
|
void sparMat::setIA(int n, int ID)
|
|
{
|
|
if (n > dim) {
|
|
cout << "ERROR in Building the sparse matrix " << endl;
|
|
exit(1);
|
|
}
|
|
rowA[n] = ID;
|
|
}
|
|
|
|
void sparMat::setJA(int n, int ID)
|
|
{
|
|
if (n >= NZ) {
|
|
cout << "ERROR " << n << " " << NZ << endl;
|
|
exit(1);
|
|
}
|
|
colA[n] = ID;
|
|
}
|
|
|
|
void sparMat::addEntry(int row, int col, fp_t val)
|
|
{
|
|
int nn;
|
|
|
|
if (type != NON_SYMMETRIC) {
|
|
if (col < row) return; // only store the upper half
|
|
nn = rowA[row];
|
|
while (colA[nn] != col && nn != rowA[row + 1]) nn ++;
|
|
if (nn == rowA[row + 1]) {
|
|
cout << "Error in Add Entry to Sparse Matrix " << endl;
|
|
cout << row << " " << col << endl;
|
|
exit(1);
|
|
}
|
|
|
|
entry[0][nn] = entry[0][nn] + val;
|
|
} else {
|
|
nn = rowA[row];
|
|
while (colA[nn] != col && nn != rowA[row + 1]) nn ++;
|
|
|
|
if (nn == rowA[row + 1]) {
|
|
cout << "Error in Add Entry to Sparse Matrix " << endl;
|
|
cout << row << " " << col << endl;
|
|
exit(1);
|
|
}
|
|
entry[0][nn] = entry[0][nn] + val; // Original version of the code
|
|
// entry[0][nn] = val; // Modified version of the code
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntry(int rDim, int cDim,
|
|
int rowMap[], int colMap[],
|
|
fp_t **entry)
|
|
{
|
|
int i, j;
|
|
for (i = 0; i < rDim; i++) {
|
|
int row = rowMap[i];
|
|
if (row < 0) continue;
|
|
|
|
for (j = 0; j < cDim; j++) {
|
|
int col = colMap[j];
|
|
if (col < 0) continue;
|
|
|
|
addEntry(row, col, entry[i][j]);
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntryNonSymm(int rDim, int cDim,
|
|
int rowMap[], int colMap[], denseMat<T> &D)
|
|
{
|
|
int i, j;
|
|
for (i = 0; i < rDim; i++) {
|
|
int row = rowMap[i];
|
|
if (row < 0) continue;
|
|
|
|
for (j = 0; j < cDim; j++) {
|
|
int col = colMap[j];
|
|
if (col < 0) continue;
|
|
fp_t check = D.getEntry(i, j);
|
|
addEntry(row, col, D.getEntry(i, j));
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntry(int nn, int mm, int *rowMap, int *colMap, denseMat<T> &D)
|
|
{
|
|
int i, j, row, col;
|
|
fp_t val;
|
|
|
|
for (i = 0; i < nn; i ++) {
|
|
row = rowMap[i];
|
|
if (row < 0) continue;
|
|
|
|
val = D.getEntry(i, i);
|
|
|
|
for (j = 0; j < mm; j ++) {
|
|
col = colMap[j];
|
|
if (col < 0) continue;
|
|
|
|
val = D.getEntry(i, j);
|
|
addEntry(row, col, val);
|
|
if (row != col) addEntry(col, row, val);
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntry(int nn, int map[], denseMat<T>& D)
|
|
{
|
|
int i, j, row, col;
|
|
fp_t val;
|
|
|
|
for (i = 0; i < nn; i ++) {
|
|
row = map[i];
|
|
if (row < 0) continue;
|
|
|
|
for (j = 0; j < nn; j ++) {
|
|
col = map[j];
|
|
if (col < 0) continue;
|
|
|
|
val = D.getEntry(i, j);
|
|
addEntry(row, col, val);
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntrySymm(int mm, int nn, int rowMap[], int colMap[],
|
|
const denseMat<T> &D)
|
|
{
|
|
int i, j, row, col;
|
|
fp_t val;
|
|
|
|
for (i = 0; i < mm; i ++) {
|
|
row = rowMap[i];
|
|
if (row < 0) continue;
|
|
|
|
for (j = 0; j < nn; j ++) {
|
|
col = colMap[j];
|
|
if (col < 0) continue;
|
|
|
|
val = D.getEntry(i, j);
|
|
if (row != col) val *= 0.5;
|
|
|
|
if (row <= col) addEntry(row, col, val);
|
|
if (row > col) addEntry(col, row, val);
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::blockAddEntryAnti(int mm, int nn, int rowMap[], int colMap[],
|
|
const denseMat<T> &D)
|
|
{
|
|
int i, j, row, col;
|
|
fp_t val;
|
|
|
|
for (i = 0; i < mm; i ++) {
|
|
row = rowMap[i];
|
|
if (row < 0) continue;
|
|
|
|
for (j = 0; j < nn; j ++) {
|
|
col = colMap[j];
|
|
if (col < 0) continue;
|
|
|
|
if (row == col) continue;
|
|
|
|
val = D.getEntry(i, j) * 0.5;
|
|
|
|
if (row < col) addEntry(row, col, val);
|
|
else {
|
|
val *= -1.;
|
|
addEntry(col, row, val);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void sparMat::diagonalScale(fp_t *d1, fp_t *d2)
|
|
{
|
|
int r;
|
|
int id;
|
|
for (r = 0; r < dim; r++) {
|
|
for (id = rowA[r]; id < rowA[r+1]; id++)
|
|
entry[0][id] *= d1[r] * d2[colA[id]];
|
|
}
|
|
}
|
|
|
|
sparMat* sparMat::symmToFull()
|
|
{
|
|
int i, j, k;
|
|
|
|
sparMat* fullMat = new sparMat;
|
|
fullMat->dim = dim;
|
|
fullMat->colDim = dim;
|
|
fullMat->NZ = NZ*2 - dim;
|
|
fullMat->NUMSEG = 1;
|
|
fullMat->SEGSIZE = NZ/NUMSEG;
|
|
fullMat->rowA = new int[dim+1];
|
|
fullMat->colA = new int[fullMat->NZ];
|
|
fullMat->entry = new fp_t*[NUMSEG];
|
|
fullMat->entry[0] = new fp_t[fullMat->NZ];
|
|
|
|
for (i = 0; i < fullMat->NZ; i++) fullMat->colA[i] = -100;
|
|
for (i = 0; i < dim+1; i++) fullMat->rowA[i] = 0;
|
|
|
|
// generate the fullMat->rowA
|
|
for (i = 0; i < dim; i++) {
|
|
for(j = rowA[i]; j < rowA[i+1]; j++) {
|
|
fullMat->rowA[i]++;
|
|
if (i < colA[j]) fullMat->rowA[colA[j]]++;
|
|
}
|
|
}
|
|
fullMat->rowA[dim] = fullMat->NZ;
|
|
for (i = dim-1; i >= 0; i--)
|
|
fullMat->rowA[i] = fullMat->rowA[i+1] - fullMat->rowA[i];
|
|
|
|
// generate the fullMat->colA
|
|
for (i = 0; i < dim; i++) {
|
|
int p_wk = fullMat->rowA[i];
|
|
|
|
while (fullMat->colA[p_wk]!=-100 && p_wk<fullMat->rowA[i+1]) {
|
|
if (p_wk == fullMat->rowA[i+1]) return 0; //error, estimation error;
|
|
p_wk++;
|
|
}
|
|
|
|
for (k = rowA[i]; k < rowA[i+1]; k++) {
|
|
j = colA[k];
|
|
|
|
fullMat->colA[p_wk] = colA[k];
|
|
fullMat->entry[0][p_wk] = entry[0][k];
|
|
p_wk++;
|
|
|
|
if (i < j) {
|
|
int pwk2 = fullMat->rowA[j];
|
|
while ((fullMat->colA[pwk2] != -100) && (pwk2 <= fullMat->rowA[j+1])) {
|
|
if (pwk2 == fullMat->rowA[j+1]) return 0; //error, estimation error;
|
|
pwk2++;
|
|
}
|
|
fullMat->colA[pwk2] = i;
|
|
fullMat->entry[0][pwk2] = entry[0][k];
|
|
}
|
|
}
|
|
}
|
|
return fullMat;
|
|
}
|
|
|
|
void sparMat::print()
|
|
{
|
|
int i, j;
|
|
|
|
cout << "dim=" << dim << endl;
|
|
cout << "colDim=" << colDim << endl;
|
|
cout << "NZ=" << NZ << endl;
|
|
cout << "SEGSIZE=" << SEGSIZE << endl;
|
|
cout << "NUMSEG=" << NUMSEG << endl;
|
|
|
|
// rowA
|
|
cout << endl << "rowA" << endl;
|
|
for (j = 0; j < dim + 1; j++)
|
|
cout << "rowA[" << j << "]=" << rowA[j] << endl;
|
|
|
|
// colA
|
|
cout << endl << "colA" << endl;;
|
|
for (j = 0; j < NZ; j++)
|
|
cout << "colA[" << j << "]=" << colA[j] << endl;
|
|
|
|
// A
|
|
for (i = 0; i < NUMSEG; i++) {
|
|
cout << endl << "A (" << i << ")" << endl;
|
|
for (j = 0; j < NZ; j++)
|
|
cout << "A[" << j << "]=" << entry[i][j] << endl;
|
|
}
|
|
|
|
// REM = NZ - NUMSEG * NZ;
|
|
// for (j = 0; j < REM; j++)
|
|
// cout << entry[NUMSEG].real() << ", " << entry[NUMSEG].imag() << endl;
|
|
|
|
}
|
|
|
|
void sparMat::print(char *fname)
|
|
{
|
|
FILE *foo;
|
|
char matName[STRLEN], rowName[STRLEN], colName[STRLEN];
|
|
|
|
// set up file names for the matrix output
|
|
memset(matName, 0, STRLEN*sizeof(char));
|
|
memset(rowName, 0, STRLEN*sizeof(char));
|
|
memset(colName, 0, STRLEN*sizeof(char));
|
|
sprintf(matName, "%s.matrix", fname);
|
|
sprintf(rowName, "%s.row", fname);
|
|
sprintf(colName, "%s.col", fname);
|
|
|
|
// col info. in binary format
|
|
if ((foo = fopen(colName, "wb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s !\n", colName);
|
|
exit(1);
|
|
}
|
|
fwrite(colA, sizeof(int), NZ, foo);
|
|
fclose(foo);
|
|
|
|
// matrix in binary format
|
|
if ((foo = fopen(matName, "wb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s ! \n", matName);
|
|
exit(1);
|
|
}
|
|
fwrite(entry[0], sizeof(fp_t), NZ, foo);
|
|
fclose(foo);
|
|
|
|
// row info. in binary format
|
|
if ((foo = fopen(rowName, "wb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s !\n", rowName);
|
|
exit(1);
|
|
}
|
|
fwrite(rowA, sizeof(int), dim + 1, foo);
|
|
fclose(foo);
|
|
}
|
|
|
|
void sparMat::printIA(char *fname)
|
|
{
|
|
FILE *foo;
|
|
char rowName[STRLEN];
|
|
memset(rowName, 0, STRLEN*sizeof(char));
|
|
sprintf(rowName, "%s.row", fname);
|
|
|
|
// row info. in binary format
|
|
if ((foo = fopen(rowName, "wb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s to Write !\n", rowName);
|
|
exit(1);
|
|
}
|
|
fwrite(rowA, sizeof(int), dim + 1, foo);
|
|
fclose(foo);
|
|
}
|
|
|
|
void sparMat::printJA(char *fname)
|
|
{
|
|
FILE *foo;
|
|
char colName[STRLEN];
|
|
memset(colName, 0, STRLEN*sizeof(char));
|
|
sprintf(colName, "%s.col", fname);
|
|
|
|
// col info. in binary format
|
|
if ((foo = fopen(colName, "wb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s !\n", colName);
|
|
exit(1);
|
|
}
|
|
fwrite(colA, sizeof(int), NZ, foo);
|
|
fclose(foo);
|
|
}
|
|
|
|
void sparMat::preRead(char *fname)
|
|
{
|
|
FILE *foo;
|
|
char rowName[STRLEN], colName[STRLEN];
|
|
|
|
// set up file names for the matrix output
|
|
memset(rowName, 0, STRLEN*sizeof(char));
|
|
memset(colName, 0, STRLEN*sizeof(char));
|
|
sprintf(rowName, "%s.row", fname);
|
|
sprintf(colName, "%s.col", fname);
|
|
|
|
// row info. in binary format
|
|
if (rowA != 0) delete [] rowA;
|
|
rowA = new int[dim + 1];
|
|
if ((foo = fopen(rowName, "rb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s to Read!\n", rowName);
|
|
exit(1);
|
|
}
|
|
fread(rowA, sizeof(int), dim + 1, foo);
|
|
fclose(foo);
|
|
NZ = rowA[dim];
|
|
|
|
// col info. in binary format
|
|
if (colA != 0) delete [] colA;
|
|
colA = new int[NZ];
|
|
if ((foo = fopen(colName, "rb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s to Read !\n", colName);
|
|
exit(1);
|
|
}
|
|
fread(colA, sizeof(int), NZ, foo);
|
|
fclose(foo);
|
|
|
|
// matrix entry
|
|
entry = new fp_t*[NUMSEG];
|
|
entry[0] = new fp_t[NZ];
|
|
}
|
|
|
|
void sparMat::read(char *fname)
|
|
{
|
|
FILE *foo;
|
|
char matName[STRLEN], rowName[STRLEN], colName[STRLEN];
|
|
|
|
// set up file names for the matrix output
|
|
memset(matName, 0, STRLEN*sizeof(char));
|
|
memset(rowName, 0, STRLEN*sizeof(char));
|
|
memset(colName, 0, STRLEN*sizeof(char));
|
|
sprintf(matName, "%s.matrix", fname);
|
|
sprintf(rowName, "%s.row", fname);
|
|
sprintf(colName, "%s.col", fname);
|
|
|
|
if (rowA == 0) rowA = new int[dim+1];
|
|
// row info. in binary format
|
|
if ((foo = fopen(rowName, "rb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s !\n", rowName);
|
|
exit(1);
|
|
}
|
|
fread(rowA, sizeof(int), dim + 1, foo);
|
|
fclose(foo);
|
|
|
|
NZ = rowA[dim];
|
|
|
|
if (colA == 0) colA = new int[NZ];
|
|
// col info. in binary format
|
|
if ((foo = fopen(colName, "rb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s !\n", colName);
|
|
exit(1);
|
|
}
|
|
fread(colA, sizeof(int), NZ, foo);
|
|
fclose(foo);
|
|
|
|
if (entry == 0) {
|
|
entry = new fp_t*[NUMSEG];
|
|
entry[0] = new fp_t[NZ];
|
|
}
|
|
// matrix in binary format
|
|
if ((foo = fopen(matName, "rb")) == NULL) {
|
|
fprintf(stderr, "Couldn't open file %s ! \n", matName);
|
|
exit(1);
|
|
}
|
|
fread(entry[0], sizeof(fp_t), NZ, foo);
|
|
fclose(foo);
|
|
}
|
|
|
|
fp_t sparMat::memUsage()
|
|
{
|
|
int mem = (dim+1) * sizeof(int) + NZ * sizeof(int) + NZ * sizeof(fp_t);
|
|
return static_cast<fp_t>(mem) / 1024. / 1024.;
|
|
}
|
|
|
|
// //////////////////////////////////////////////
|
|
// // free functions
|
|
//
|
|
// // sparse matrix matrix vector multipication
|
|
// Vector MxV(sparMat *const A, Vector &x)
|
|
// {
|
|
// int i, cc, id;
|
|
// int dim = A->getRowDim();
|
|
// fp_t cval;
|
|
//
|
|
//
|
|
// if (A->getColDim() != static_cast<int>(x.size())) {
|
|
// cout << " can not match Matrix-vector dimensions in dd_sparmat::MxV"
|
|
// << endl;
|
|
// cout << "A->getColDim=" << A->getColDim() << " x.size=" << x.size() <<endl;
|
|
// exit(0);
|
|
// }
|
|
//
|
|
// Vector y((unsigned int)A->getRowDim());
|
|
// for (i = 0; i < dim; i++) {
|
|
// for (id = A->rowA[i]; id < A->rowA[i + 1]; id++) {
|
|
// cc = A->colA[id];
|
|
// cval = A->getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
//
|
|
// // // debug
|
|
// // printf("i=%3d x[%4d]=(%8.4f, %8.4f) A[%4d]=(%8.4f, %8.4f)\n",
|
|
// // i, cc, x[cc].real(), x[cc].imag(), id, A->getEntry(id).real(),
|
|
// // A->getEntry(id).imag());
|
|
// }
|
|
// }
|
|
//
|
|
// return y;
|
|
// }
|
|
//
|
|
// Vector MxV(sparMat *const A, Vector &x, int *const xMap, int *const yMap)
|
|
// {
|
|
// int i, cc, id;
|
|
// int dimR = A->getRowDim();
|
|
// fp_t cval;
|
|
//
|
|
// if (A->getColDim() != static_cast<int>(x.size())) {
|
|
// cout << "can not match Matrix-vector dimensions in dd_sparmat::MxV" <<endl;
|
|
// cout << "A->getColDim=" << A->getColDim() << " x.size=" << x.size() <<endl;
|
|
// exit(0);
|
|
// }
|
|
//
|
|
// Vector y((unsigned int)dimR);
|
|
// for (i = 0; i < dimR; i++) {
|
|
// for (id = A->rowA[i]; id < A->rowA[i+1]; id++) {
|
|
// cc = A->colA[id];
|
|
// cval = A->getEntry(id) * x[ xMap[cc] ];
|
|
// y[ yMap[i+dimR] ] += cval;
|
|
// }
|
|
// }
|
|
//
|
|
// return y;
|
|
// }
|
|
//
|
|
// // sparse matrix and transpose matrix vector multipication
|
|
// Vector MxV(sparMat *const A, Vector &x, bool transpose)
|
|
// {
|
|
// Vector y((unsigned int)A->getDim());
|
|
//
|
|
// if (A->isSymmetric() == NON_SYMMETRIC)
|
|
// nonSymmMxV(A, x, y, transpose);
|
|
// else
|
|
// symmMxV(A, x, y);
|
|
//
|
|
// return y;
|
|
// }
|
|
//
|
|
// // symmetric sparse matrix and transpose matrix vector multiplication
|
|
// void symmMxV(sparMat *const A, Vector &x, Vector &y)
|
|
// {
|
|
// if (A->getColDim() != static_cast<int>(x.size())) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A->getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A->isSymmetric() == NON_SYMMETRIC) {
|
|
// cout << "matrix is non-symmtric symmMxV can not be used in this case" << endl;
|
|
// exit(0);
|
|
// }
|
|
// int i, cc, id;
|
|
// int dim = A->getDim();
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// id = A->rowA[i];
|
|
// cval=A->getEntry(id) * x[i] ;
|
|
// y[i] += cval;
|
|
// for (id = A->rowA[i] + 1; id < A->rowA[i + 1]; id++) {
|
|
// cc = A->colA[id];
|
|
// if (cc < dim) {
|
|
// cval = A->getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// cval = A->getEntry(id) * x[i];
|
|
// y[cc] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// }
|
|
//
|
|
// // non symetric sparse matrix and transpose matrix vector multipication
|
|
void sparMat::NonSymmMxV(fp_t *x, fp_t *y)
|
|
{
|
|
// int i, cc, id;
|
|
|
|
// first reset y
|
|
#pragma omp parallel for schedule(static)
|
|
for (int i = 0; i < dim ; i ++) y[i] = 0.0;
|
|
|
|
// fp_t rval;
|
|
#pragma omp parallel for schedule(static)
|
|
for (int i = 0; i < dim; i++) {
|
|
for (int id = rowA[i]; id < rowA[i + 1]; id ++) {
|
|
int cc = colA[id];
|
|
fp_t rval = getEntry(id) * x[cc];
|
|
y[i] += rval;
|
|
// if(abs(y[i])>1.0E-16) cout<<y[i]<<endl;
|
|
}
|
|
// cout << "y[" << i << "] = " << y[i] << endl;
|
|
}
|
|
}
|
|
//
|
|
|
|
|
|
// void nonSymmMxV(sparMat *const A, Vector &x, Vector &y, bool transpose)
|
|
// {
|
|
// if (!transpose) {
|
|
// if (A->getColDim() != static_cast<int>(x.size())) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A col dim=" << A->getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A->getRowDim() != static_cast<int>(y.size())) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "y size=" << y.size() << " A row dim=" << A->getDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// } else {
|
|
// if (A->getRowDim() != static_cast<int>(x.size())) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " transpose(A) col dim=" << A->getDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A->getColDim() != static_cast<int>(y.size())) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "y size=" << y.size() << " transpose(A) row dim=" << A->getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// }
|
|
//
|
|
// if (A->isSymmetric() != NON_SYMMETRIC) {
|
|
// cout << "matrix isn't non-symmetric nonSymmMxV cannot be used in this case"
|
|
// << endl;
|
|
// exit(0);
|
|
// }
|
|
//
|
|
// int i, cc, id;
|
|
// int dim = A->getDim();
|
|
// fp_t cval;
|
|
//
|
|
// if (transpose) {
|
|
// for (i = 0; i < dim; i++) {
|
|
// for (id = A->rowA[i]; id < A->rowA[i + 1]; id++) {
|
|
// cc = A->colA[id];
|
|
// cval = A->getEntry(id) * x[i];
|
|
// y[cc] += cval;
|
|
// }
|
|
// }
|
|
// } else {
|
|
// for (i = 0; i < dim; i++) {
|
|
// for (id = A->rowA[i]; id < A->rowA[i + 1]; id++) {
|
|
// cc = A->colA[id];
|
|
// cval = A->getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // symmetric matrix-vector multiplication
|
|
// void symm_MxV(const sparMat &A, Vector &x, Vector &y)
|
|
// {
|
|
// if (static_cast<int>(x.size()) != A.getColDim()) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A.getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A.type != SYMMETRIC) {
|
|
// cout << "The matrix is not symmtric. symm_MxV can not be used in this case"
|
|
// << endl;
|
|
// exit(0);
|
|
// }
|
|
// int i, cc, id;
|
|
// int dim = A.getDim();
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// id = A.getRowMap(i);
|
|
// cval=A.getEntry(id) * x[i] ;
|
|
// y[i] += cval;
|
|
// for (id = A.getRowMap(i) + 1; id < A.getRowMap(i+1); id ++) {
|
|
// cc = A.getColMap(id);
|
|
// if (cc < dim) {
|
|
// cval = A.getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// cval = A.getEntry(id) * x[i];
|
|
// y[cc] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // non-symmetric sparse matrix matrix-vector multipication
|
|
// void nonsymm_MxV(sparMat const &A, Vector &x, Vector &y)
|
|
// {
|
|
// if (x.size() != static_cast<unsigned int>(A.getColDim()) ) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A.getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
//
|
|
// if (A.type == SYMMETRIC ) {
|
|
// cout << "matrix is symmtric nonSymmMxV can not be used in this case"
|
|
// << endl;
|
|
// exit(0);
|
|
// }
|
|
//
|
|
// int i, cc, id;
|
|
// int dim = A.getDim();
|
|
//
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// for (id = A.rowA[i]; id < A.rowA[i + 1]; id ++) {
|
|
// cc = A.colA[id];
|
|
// cval = A.getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // anti-symmetric matrix-vector multiplication
|
|
// void antisymm_MxV(const sparMat &A, Vector &x, Vector &y)
|
|
// {
|
|
// if (static_cast<int>(x.size()) != A.getColDim()) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A.getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A.type != ANTI_SYMMETRIC) {
|
|
// cout << "The matrix is not anti-symmtric. antisymm_MxV can not be used in this case" << endl;
|
|
// exit(0);
|
|
// }
|
|
// int i, cc, id;
|
|
// int dim = A.getDim();
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// id = A.getRowMap(i);
|
|
// cval=A.getEntry(id) * x[i] ;
|
|
// y[i] += cval;
|
|
// for (id = A.getRowMap(i) + 1; id < A.getRowMap(i+1); id ++) {
|
|
// cc = A.getColMap(id);
|
|
// if (cc < dim) {
|
|
// cval = A.getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// cval = A.getEntry(id) * x[i];
|
|
// y[cc] -= cval;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // hermitian matrix-vector multiplication
|
|
// void hermitian_MxV(const sparMat &A, Vector &x, Vector &y)
|
|
// {
|
|
// if (static_cast<int>(x.size()) != A.getColDim()) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A.getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A.type != HERMITIAN) {
|
|
// cout << "The matrix is not hermitian. antisymm_MxV can not be used in this case" << endl;
|
|
// exit(0);
|
|
// }
|
|
// int i, cc, id;
|
|
// int dim = A.getDim();
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// id = A.getRowMap(i);
|
|
// cval=A.getEntry(id) * x[i] ;
|
|
// y[i] += cval;
|
|
// for (id = A.getRowMap(i) + 1; id < A.getRowMap(i+1); id ++) {
|
|
// cc = A.getColMap(id);
|
|
// if (cc < dim) {
|
|
// cval = A.getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// cval = conj(A.getEntry(id)) * x[i];
|
|
// y[cc] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // skew-hermitian matrix-vector multiplication
|
|
// void skew_herm_MxV(const sparMat &A, Vector &x, Vector &y)
|
|
// {
|
|
// if (static_cast<int>(x.size()) != A.getColDim()) {
|
|
// cout << "incompatble matrix and vector dimensions in sparse MxV" << endl;
|
|
// cout << "x size=" << x.size() << " A dim=" << A.getColDim() << endl;
|
|
// exit(0);
|
|
// }
|
|
// if (A.type != SKEW_HERMITIAN) {
|
|
// cout << "The matrix is not skew-hermitian. antisymm_MxV can not be used in this case" << endl;
|
|
// exit(0);
|
|
// }
|
|
// int i, cc, id;
|
|
// int dim = A.getDim();
|
|
// fp_t cval;
|
|
//
|
|
// for (i = 0; i < dim; i++) {
|
|
// id = A.getRowMap(i);
|
|
// cval=A.getEntry(id) * x[i] ;
|
|
// y[i] += cval;
|
|
// for (id = A.getRowMap(i) + 1; id < A.getRowMap(i+1); id ++) {
|
|
// cc = A.getColMap(id);
|
|
// if (cc < dim) {
|
|
// cval = A.getEntry(id) * x[cc];
|
|
// y[i] += cval;
|
|
// cval = -conj(A.getEntry(id)) * x[i];
|
|
// y[cc] += cval;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// Vector operator*(const sparMat &A, Vector &x)
|
|
// {
|
|
// int dim = A.getDim();
|
|
// Vector y((unsigned int)dim);
|
|
//
|
|
// MatrixType type = A.type;
|
|
//
|
|
// switch (type) {
|
|
// case SYMMETRIC:
|
|
// // cout << "symmetric matrix-vector multiplication" << endl;
|
|
// symm_MxV(A, x, y);
|
|
// break;
|
|
// case NON_SYMMETRIC:
|
|
// // cout << "non-symmetric matrix-vector multiplication" << endl;
|
|
// nonsymm_MxV(A, x, y);
|
|
// break;
|
|
// case ANTI_SYMMETRIC:
|
|
// // cout << "anti-symmetric matrix-vector multiplication" << endl;
|
|
// antisymm_MxV(A, x, y);
|
|
// break;
|
|
// case HERMITIAN:
|
|
// // cout << "hermitian matrix-vector multiplication" << endl;
|
|
// hermitian_MxV(A, x, y);
|
|
// break;
|
|
// case SKEW_HERMITIAN:
|
|
// // cout << "skew-hermitian matrix-vector multiplication" << endl;
|
|
// skew_herm_MxV(A, x, y);
|
|
// break;
|
|
// }
|
|
//
|
|
// return y;
|
|
// }
|
|
|
|
void sparMat::st_to_hb(int matdim,int rowdim, int nst,int ist[],int jst[],fp_t vst[])
|
|
{
|
|
int row, col, cnt(0),pre(0),CNT(0);
|
|
fp_t val;
|
|
|
|
NZ = nst;
|
|
setDim(matdim);
|
|
setColDim(rowdim);//not important just write to header
|
|
|
|
setMatType(NON_SYMMETRIC);
|
|
SEGSIZE = 1;
|
|
|
|
entry = new fp_t*[NUMSEG];
|
|
colA = new int[NZ];
|
|
rowA = ( int * ) malloc ( ( matdim + 1 ) * sizeof ( int ) );
|
|
|
|
int ncc= st_to_cc_size ( nst, ist, jst );
|
|
printf ( "\n" );
|
|
printf ( " Number of CC values = %d\n", ncc );
|
|
printf ( " Matrix dimension is = %d\n", matdim );
|
|
/*
|
|
Create the CC indices.
|
|
*/
|
|
|
|
st_to_cc_index ( nst, ist, jst, ncc, matdim, colA, rowA );
|
|
/*
|
|
Create the CC values.
|
|
*/
|
|
entry[0] = st_to_cc_values ( nst, ist, jst, vst, ncc, matdim, colA, rowA );
|
|
}
|
|
|
|
void sparMat::hb_file_write(char const * fname)
|
|
{
|
|
string output_file = fname;
|
|
ofstream output;
|
|
int totcrd,ptrcrd,indcrd,valcrd,rhscrd,nrow,ncol,nnzero,neltvl,nrhs,nrhsix;
|
|
char* MXTYPE;
|
|
char title[73] = "1Real unsymmetric assembled matrix based on DGTD";
|
|
valcrd = ceil(NZ/10);
|
|
char key[9] = "RUA_32";
|
|
char mxtype[4] = "RUA";
|
|
char ptrfmt[17] = "(13I8)";
|
|
char indfmt[17] = "(13I8)";
|
|
char valfmt[21] = "(10E16.8)";//(10E16.8)
|
|
char rhsfmt[21] = "(10E16.8)";
|
|
|
|
cout << " HB_FILE_WRITE writes an HB file.\n";
|
|
cout << "\n";
|
|
cout << " Writing the file '" << output_file << "'.\n";
|
|
|
|
neltvl = nrhs = nrhsix = 0;
|
|
nnzero = NZ;
|
|
ncol = dim;//ncol is real number of "row"
|
|
nrow = colDim;//nrow just been writen to file header
|
|
ptrcrd = ceil((dim+1)/13);
|
|
indcrd = ceil(nnzero/13);
|
|
valcrd = ceil(nnzero/10);//10
|
|
rhscrd = 1;//special
|
|
totcrd = ptrcrd+indcrd+valcrd;
|
|
|
|
output.open ( output_file.c_str ( ) );
|
|
if(!output){
|
|
cout << "\n TEST10 - Warning!\n Error opening the file.\n";
|
|
return;
|
|
}
|
|
hb_file_write_no_rhs(output, title, key, totcrd, ptrcrd, indcrd,
|
|
valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt,
|
|
valfmt, rhsfmt, NULL, nrhs, nrhsix, rowA, colA, entry[0],
|
|
NULL, NULL, NULL, NULL, NULL , NULL );
|
|
|
|
output.close ( );
|
|
|
|
return;
|
|
}
|
|
|
|
|