Loading files
This file is pinned to the selected release. Download the original to open it in your design tool.
// GLKSparseMatrix.cpp: implementation of the GLKSparseMatrix class.
//
//////////////////////////////////////////////////////////////////////
#define _CRT_SECURE_NO_DEPRECATE
#include <stdio.h>
#include "GLKHeap.h"
#include "GLKObList.h"
#include "GLKSparseMatrix.h"
/*//-------------------------------------------------------------------
// Code in the following lines are for detecting memory leaks
#ifdef _DEBUG
#include <afx.h>
#define new DEBUG_NEW
#endif*/
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
GLKSparseMatrix::GLKSparseMatrix()
{
m_nCol=0; m_nRow=0;
}
GLKSparseMatrix::~GLKSparseMatrix()
{
_clearAll();
}
//////////////////////////////////////////////////////////////////////
// Implementation
//////////////////////////////////////////////////////////////////////
void GLKSparseMatrix::Initialization(unsigned long nRow, unsigned long nCol)
{
unsigned long i;
_clearAll();
// if (nCol>nRow) return;
if (nRow>0) {
m_nRow=nRow; m_nCol=nCol;
m_rowRoots=(GLKSparseMatrixNode**)new long[m_nRow];
m_rowElementNum=new unsigned long[m_nRow];
for(i=0;i<m_nRow;i++) {m_rowRoots[i]=NULL; m_rowElementNum[i]=0;}
m_permutations=new unsigned long[m_nRow];
for(i=0;i<m_nRow;i++) m_permutations[i]=i;
}
}
void GLKSparseMatrix::InputElement(unsigned long nRowIndex, unsigned long nColIndex, double data)
{
GLKSparseMatrixNode *newNode=new GLKSparseMatrixNode;
newNode->nColIndex=nColIndex; newNode->data=data;
newNode->pNext=m_rowRoots[nRowIndex];
m_rowRoots[nRowIndex]=newNode;
m_rowElementNum[nRowIndex]++;
}
void GLKSparseMatrix::AddElement(unsigned long nRowIndex, unsigned long nColIndex, double delta)
{
if(m_rowRoots[nRowIndex]){
GLKSparseMatrixNode* currentNode = m_rowRoots[nRowIndex];
for(;currentNode!=NULL;){
GLKSparseMatrixNode* nextNode=currentNode->pNext;
if(currentNode->nColIndex==nColIndex){
currentNode->data+=delta;
return;
}
currentNode=nextNode;
}
}
InputElement(nRowIndex,nColIndex,delta);
}
void GLKSparseMatrix::EliminateRow(unsigned long nRowIndex)
{
GLKSparseMatrixNode *currentNode=m_rowRoots[nRowIndex];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
delete currentNode;
currentNode=nextNode;
}
m_rowRoots[nRowIndex]=NULL;
m_rowElementNum[nRowIndex]=0;
}
unsigned long GLKSparseMatrix::GetNonzeroElementNumber()
{
unsigned long i,nnzNum=0;
for(i=0;i<m_nRow;i++) nnzNum+=m_rowElementNum[i];
return nnzNum;
}
void GLKSparseMatrix::InputComplete(bool bSort)
{
unsigned long i,j,num;
//-------------------------------------------------------------------------
// resort the elements in each row
for(i=0;i<m_nRow;i++) {
GLKSparseMatrixNode *currentNode=m_rowRoots[i];
num=m_rowElementNum[i]; if (num==0) continue;
GLKSparseMatrixNode **nodeArray=(GLKSparseMatrixNode **)new long[num];
if (bSort) {
currentNode=m_rowRoots[i];
j=0;
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
nodeArray[j]=currentNode;
currentNode=nextNode; j++;
}
_quickSort(nodeArray,0,num-1,num);
}
else {
currentNode=m_rowRoots[i];
j=0;
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
nodeArray[num-1-j]=currentNode;
currentNode=nextNode; j++;
}
}
for(j=0;j<num-1;j++) nodeArray[j]->pNext=nodeArray[j+1];
nodeArray[num-1]->pNext=NULL; m_rowRoots[i]=nodeArray[0];
// for(j=0;j<num;j++) printf("%ld (%lf) ",nodeArray[j]->nColIndex,nodeArray[j]->data); printf("\n");
delete [](GLKSparseMatrixNode **)nodeArray;
}
}
void GLKSparseMatrix::InputPermutationTable(unsigned long *permutations)
{
unsigned long i;
for(i=0;i<m_nRow;i++) m_permutations[i]=permutations[i];
}
void GLKSparseMatrix::OutputPermutationTable(unsigned long *permutations)
{
unsigned long i;
for(i=0;i<m_nRow;i++) permutations[i]=m_permutations[i];
}
void GLKSparseMatrix::ReOrderingByCuthillMcKeeAlgorithm(bool bSymmetricMatrix)
{
if (m_nCol!=m_nRow) return; // the algorithm can only be applied to the square matrix
bool *bInserted; unsigned long *rSet; unsigned long rSetNum=0;
unsigned long i;
bInserted=new bool[m_nRow]; rSet=new unsigned long[m_nRow];
for(i=0;i<m_nRow;i++) bInserted[i]=false;
GLKSparseMatrixNode** m_colRoots;
GLKSparseMatrixNode** m_colTails;
if (!bSymmetricMatrix) {
m_colRoots=(GLKSparseMatrixNode**)new long[m_nCol];
m_colTails=(GLKSparseMatrixNode**)new long[m_nCol];
for(i=0;i<m_nCol;i++) {m_colTails[i]=NULL;m_colRoots[i]=NULL;}
for(i=0;i<m_nRow;i++) {
GLKSparseMatrixNode *currentNode=m_rowRoots[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
GLKSparseMatrixNode *newNode=new GLKSparseMatrixNode;
newNode->nColIndex=i;
unsigned long colIndex=currentNode->nColIndex;
if (m_colTails[colIndex]) {
m_colTails[colIndex]->pNext=newNode;
m_colTails[colIndex]=newNode;
}
else {
m_colRoots[colIndex]=newNode;
m_colTails[colIndex]=newNode;
}
currentNode=nextNode;
}
}
delete [](GLKSparseMatrixNode**)m_colTails;
}
//--------------------------------------------------------------------------------------
// Step 1: find the peripheral vertex and insert it into the rSet
unsigned long peripheralVertexIndex;
unsigned long degree=m_nRow+1;
for(i=0;i<m_nRow;i++) {
if (m_rowElementNum[i]<degree) { degree=m_rowElementNum[i]; peripheralVertexIndex=i; }
}
// peripheralVertexIndex=_pseudoPeripheralNodeSearch();
rSet[rSetNum++]=peripheralVertexIndex; bInserted[peripheralVertexIndex]=true;
//--------------------------------------------------------------------------------------
// Step 2: incrementally add more nodes into the rSet
for(i=0;rSetNum<m_nRow;i++) {
GLKSparseMatrixNode *currentNode=m_rowRoots[rSet[i]];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
if (!(bInserted[currentNode->nColIndex])) {
rSet[rSetNum++]=currentNode->nColIndex;
bInserted[currentNode->nColIndex]=true;
}
currentNode=nextNode;
}
if (i==rSetNum-1) { // insert another peripheral vertex
unsigned long peripheralVertexIndex;
unsigned long degree=m_nRow+1;
for(unsigned j=0;j<m_nRow;j++) {
if (bInserted[j]) continue;
if (m_rowElementNum[j]<degree) { degree=m_rowElementNum[j]; peripheralVertexIndex=j; }
}
rSet[rSetNum++]=peripheralVertexIndex; bInserted[peripheralVertexIndex]=true;
}
if (bSymmetricMatrix) continue;
currentNode=m_colRoots[rSet[i]];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
if (!(bInserted[currentNode->nColIndex])) {
rSet[rSetNum++]=currentNode->nColIndex;
bInserted[currentNode->nColIndex]=true;
}
currentNode=nextNode;
}
}
//--------------------------------------------------------------------------------------
// Step 3: build the permulation table
for(i=0;i<rSetNum;i++) m_permutations[rSet[i]]=i;
// for(i=0;i<rSetNum;i++) m_permutations[rSet[rSetNum-1-i]]=i;
if (!bSymmetricMatrix) {
for(i=0;i<m_nCol;i++) {
GLKSparseMatrixNode *currentNode=m_colRoots[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
delete currentNode;
currentNode=nextNode;
}
}
delete [](GLKSparseMatrixNode**)m_colRoots;
}
delete []bInserted; delete []rSet;
}
int GLKSparseMatrix::_pseudoPeripheralNodeSearch(GLKSparseMatrixNode** colRoots)
{
int levelNum,newRootIndex,levelNum2,newRootIndex2;
_levelStructuring((int)m_nRow/2,levelNum,newRootIndex,colRoots);
while(true) {
_levelStructuring(newRootIndex,levelNum2,newRootIndex2,colRoots);
if (levelNum2<=levelNum) return newRootIndex;
newRootIndex=newRootIndex2; levelNum=levelNum2;
}
return 0;
}
void GLKSparseMatrix::_levelStructuring(int rootIndex, int &levelNum, int &newRootIndex,
GLKSparseMatrixNode** colRoots)
{
unsigned long i,currentLevel;
GLKArray *bottomLevelIndexArray,*currentLevelIndexArray;
bottomLevelIndexArray=new GLKArray(100,100,1);
currentLevelIndexArray=new GLKArray(100,100,1);
int* nodeLevels=new int[m_nRow];
for(i=0;i<m_nRow;i++) nodeLevels[i]=-1;
//-------------------------------------------------------------------------
// Build the level structure
currentLevelIndexArray->Add(rootIndex); nodeLevels[rootIndex]=0;
currentLevel=0;
do{
bottomLevelIndexArray->RemoveAll();
int j,num=currentLevelIndexArray->GetSize();
for(j=0;j<num;j++) {
int elementIndex=currentLevelIndexArray->GetIntAt(j);
GLKSparseMatrixNode *currentNode=m_rowRoots[elementIndex];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
if (nodeLevels[currentNode->nColIndex]<0) {
nodeLevels[currentNode->nColIndex]=currentLevel+1;
bottomLevelIndexArray->Add((int)(currentNode->nColIndex));
}
currentNode=nextNode;
}
if (!colRoots) continue;
currentNode=colRoots[elementIndex];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
if (nodeLevels[currentNode->nColIndex]<0) {
nodeLevels[currentNode->nColIndex]=currentLevel+1;
bottomLevelIndexArray->Add((int)(currentNode->nColIndex));
}
currentNode=nextNode;
}
}
if (bottomLevelIndexArray->GetSize()==0) break;
currentLevelIndexArray->RemoveAll();
num=bottomLevelIndexArray->GetSize();
for(j=0;j<num;j++) currentLevelIndexArray->Add(bottomLevelIndexArray->GetIntAt(j));
currentLevel++;
}while(true);
//-------------------------------------------------------------------------
// Find the one with the minimum degree in the bottom level
int num=currentLevelIndexArray->GetSize(); unsigned long degree=m_nRow;
for(int j=0;j<num;j++) {
int index=currentLevelIndexArray->GetIntAt(j);
if (m_rowElementNum[index]<degree) {
degree=m_rowElementNum[index]; newRootIndex=index;
}
}
levelNum=currentLevel+1;
delete bottomLevelIndexArray; delete currentLevelIndexArray;
delete []nodeLevels;
}
void GLKSparseMatrix::_quickSort(GLKSparseMatrixNode **nodeArray,
long left, long right, long max)
{
long i,j; GLKSparseMatrixNode *s;
if (left < right){
s = nodeArray[left];
i = left; j = right + 1;
while(true) {
while( (i + 1 < max) && (nodeArray[++i]->nColIndex < s->nColIndex) ) ;
while( (j -1 > -1) && (nodeArray[--j]->nColIndex > s->nColIndex) ) ;
if (i >= j) break;
GLKSparseMatrixNode * temp=nodeArray[i];
nodeArray[i]=nodeArray[j];
nodeArray[j]=temp;
}
nodeArray[left] = nodeArray[j];
nodeArray[j] = s;
_quickSort(nodeArray, left, j-1, max);
_quickSort(nodeArray, j+1, right, max);
}
}
void GLKSparseMatrix::_clearAll()
{
unsigned long i;
if ((m_nRow==0) && (m_nCol==0)) return;
for(i=0;i<m_nRow;i++) {
GLKSparseMatrixNode *currentNode=m_rowRoots[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
delete currentNode;
currentNode=nextNode;
}
}
delete [](GLKSparseMatrixNode **)m_rowRoots;
delete []m_rowElementNum;
delete []m_permutations;
m_nRow=0; m_nCol=0;
}
void GLKSparseMatrix::DEBUG_OutputPattern(char *filename)
{
FILE *fp; unsigned long i,j,fnum; double xx,yy;
fp=fopen(filename, "w");
fnum=0;
for(i=0;i<m_nRow;i++) {
GLKSparseMatrixNode *currentNode=m_rowRoots[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
j=currentNode->nColIndex;
xx=0.1*(double)(m_permutations[j])-0.1*(double)(m_nCol)*0.5;
yy=0.1*(double)(m_nRow-1-m_permutations[i])-0.1*(double)(m_nRow)*0.5;
fprintf(fp,"v %lf %lf 0.0\n",xx,yy);
currentNode=nextNode; fnum++;
}
}
fprintf(fp,"v %lf %lf -1.0\n",-0.1*(double)(m_nCol)*0.5,-0.1*(double)(m_nRow)*0.5);
fprintf(fp,"v %lf %lf -1.0\n",-0.1*(double)(m_nCol)*0.5,0.1*(double)(m_nRow)-0.1*(double)(m_nRow)*0.5);
fprintf(fp,"v %lf %lf -1.0\n",0.1*(double)m_nCol*0.5,0.1*(double)(m_nRow)-0.1*(double)(m_nRow)*0.5);
fprintf(fp,"v %lf %lf -1.0\n",0.1*(double)m_nCol*0.5,-0.1*(double)(m_nRow)*0.5);
fprintf(fp,"f %ld %ld %ld %ld\n",fnum+1,fnum+2,fnum+3,fnum+4);
fclose(fp);
}
void Mul_A_b(GLKSparseMatrix *sparseMat, double *b, double *x)
{
unsigned long i;
unsigned long nRow=sparseMat->GetRowNumber();
for(i=0;i<nRow;i++) {
x[i]=0.0;
GLKSparseMatrixNode* currentNode=(sparseMat->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
x[i]+=b[currentNode->nColIndex]*(currentNode->data);
currentNode=nextNode;
}
}
}
void Mul_At_b(GLKSparseMatrix *sparseMat, double *b, double *x)
{
unsigned long i,nRow;
GLKSparseMatrix *outputMat;
Transpose(sparseMat,outputMat);
nRow=outputMat->GetRowNumber();
for(i=0;i<nRow;i++) {
x[i]=0.0;
GLKSparseMatrixNode* currentNode=(outputMat->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
x[i]+=b[currentNode->nColIndex]*(currentNode->data);
currentNode=nextNode;
}
}
delete outputMat;
}
void Mul_At_A1(GLKSparseMatrix *sparseMat, GLKSparseMatrix *&outputMat)
{
unsigned long i, nRow, nCol, nRow1;
outputMat=new GLKSparseMatrix; nRow = nCol = sparseMat->GetColNumber();
nRow1 = sparseMat->GetRowNumber();
outputMat->Initialization(nRow,nCol);
for(i=0;i<nRow1;i++){
GLKSparseMatrixNode* currentnode1 = (sparseMat->GetRowRoots())[i];
for(;currentnode1!=NULL;){
GLKSparseMatrixNode* currentnode2 = (sparseMat->GetRowRoots())[i];
for(;currentnode2!=NULL;){
outputMat->AddElement(currentnode1->nColIndex, currentnode2->nColIndex, (currentnode1->data)*(currentnode2->data));
currentnode2 = currentnode2->pNext;
}
currentnode1 = currentnode1->pNext;
}
}
outputMat->InputComplete(false);
}
void Mul_At_A(GLKSparseMatrix *sparseMat, GLKSparseMatrix *&outputMat)
{
unsigned long i,j,nRow,nCol; double data;
GLKSparseMatrix *transposeMat;
Transpose(sparseMat,transposeMat);
outputMat=new GLKSparseMatrix; nRow=nCol=sparseMat->GetColNumber();
outputMat->Initialization(nRow,nCol);
for(i=0;i<nRow;i++) {
for(j=0;j<nCol;j++) {
GLKSparseMatrixNode* currentRowNode=(transposeMat->GetRowRoots())[i];
GLKSparseMatrixNode* currentColNode=(transposeMat->GetRowRoots())[j];
data=0.0;
for(;(currentRowNode!=NULL)&&(currentColNode!=NULL);) {
if (currentRowNode->nColIndex==currentColNode->nColIndex) {
data+=(currentRowNode->data)*(currentColNode->data);
currentRowNode=currentRowNode->pNext;
currentColNode=currentColNode->pNext;
}
else if (currentRowNode->nColIndex<currentColNode->nColIndex) {
currentRowNode=currentRowNode->pNext;
}
else {
currentColNode=currentColNode->pNext;
}
}
if (data!=0.0) outputMat->InputElement(i,j,data);
}
}
outputMat->InputComplete(false);
delete transposeMat;
}
void Mul_A_At(GLKSparseMatrix *sparseMat, GLKSparseMatrix *&transposeMat, GLKSparseMatrix *&outputMat)
{
unsigned long i, nRow, nCol, nRow1;
Transpose(sparseMat,transposeMat);
outputMat=new GLKSparseMatrix; nRow=nCol=transposeMat->GetColNumber();
nRow1 = transposeMat->GetRowNumber();
outputMat->Initialization(nRow,nCol);
for(i=0;i<nRow1;i++){
GLKSparseMatrixNode* currentnode1 = (transposeMat->GetRowRoots())[i];
for(;currentnode1!=NULL;){
GLKSparseMatrixNode* currentnode2 = (transposeMat->GetRowRoots())[i];
for(;currentnode2!=NULL;){
outputMat->AddElement(currentnode1->nColIndex, currentnode2->nColIndex, (currentnode1->data)*(currentnode2->data));
currentnode2 = currentnode2->pNext;
}
currentnode1 = currentnode1->pNext;
}
}
outputMat->InputComplete(false);
}
void CopyMat(GLKSparseMatrix *sparseMat, GLKSparseMatrix *&outputMat)
{
unsigned long nRow,nCol,i;
nRow=sparseMat->GetRowNumber(); nCol=sparseMat->GetColNumber();
outputMat=new GLKSparseMatrix; outputMat->Initialization(nRow,nCol);
for(i=0;i<nRow;i++) {
GLKSparseMatrixNode* currentNode=(sparseMat->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
outputMat->InputElement(i,currentNode->nColIndex,currentNode->data);
currentNode=nextNode;
}
}
outputMat->InputComplete(false);
}
void PrintMat(GLKSparseMatrix *sparseMat)
{
int nRowNum=sparseMat->GetRowNumber();
for(int i=0;i<nRowNum;i++) {
GLKSparseMatrixNode *currentNode=(sparseMat->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode *nextNode=currentNode->pNext;
printf("rowIndex=%d colIndex=%d value=%lf\n",i,currentNode->nColIndex,currentNode->data);
currentNode=nextNode;
}
}
}
bool Mul_A_B(GLKSparseMatrix *matA, GLKSparseMatrix *matB, GLKSparseMatrix *&outputMat)
{
unsigned long aRow,aCol,bRow,bCol,i,j; double data;
aRow=matA->GetRowNumber(); aCol=matA->GetColNumber();
bRow=matB->GetRowNumber(); bCol=matB->GetColNumber();
if (aCol!=bRow) {outputMat=NULL; return false;}
//--------------------------------------------------------------------------------
// Preparation
GLKSparseMatrixNode **colRootNodes; colRootNodes=(GLKSparseMatrixNode **)new long[bCol];
GLKSparseMatrixNode **colTailNodes; colTailNodes=(GLKSparseMatrixNode **)new long[bCol];
for(j=0;j<bCol;j++) colRootNodes[j]=colTailNodes[j]=NULL;
for(i=0;i<bRow;i++) {
GLKSparseMatrixNode* currentNode=(matB->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
j=currentNode->nColIndex;
GLKSparseMatrixNode* newNode=new GLKSparseMatrixNode;
newNode->nColIndex=i;
newNode->data=currentNode->data;
if (colTailNodes[j]) {
colTailNodes[j]->pNext=newNode; colTailNodes[j]=newNode;
}
else {
colRootNodes[j]=newNode; colTailNodes[j]=newNode;
}
currentNode=nextNode;
}
}
//--------------------------------------------------------------------------------
// Multiplication
outputMat=new GLKSparseMatrix; outputMat->Initialization(aRow,bCol);
for(i=0;i<aRow;i++) {
for(j=0;j<bCol;j++) {
GLKSparseMatrixNode* currentRowNode=(matA->GetRowRoots())[i];
GLKSparseMatrixNode* currentColNode=colRootNodes[j];
data=0.0; bool bFlag=false;
for(;(currentRowNode!=NULL)&&(currentColNode!=NULL);) {
if (currentRowNode->nColIndex==currentColNode->nColIndex) {
data+=(currentRowNode->data)*(currentColNode->data);
currentRowNode=currentRowNode->pNext;
currentColNode=currentColNode->pNext;
bFlag=true;
}
else if (currentRowNode->nColIndex<currentColNode->nColIndex) {
currentRowNode=currentRowNode->pNext;
}
else {
currentColNode=currentColNode->pNext;
}
}
if (bFlag) outputMat->InputElement(i,j,data);
}
}
outputMat->InputComplete(false);
//--------------------------------------------------------------------------------
// Free the memory
for(j=0;j<bCol;j++) {
GLKSparseMatrixNode* currentNode=colRootNodes[j];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
delete currentNode;
currentNode=nextNode;
}
}
delete [](GLKSparseMatrixNode **)colRootNodes; delete [](GLKSparseMatrixNode **)colTailNodes;
return true;
}
void Transpose(GLKSparseMatrix *sparseMat, GLKSparseMatrix *&outputMat)
{
unsigned long i,nRow;
outputMat=new GLKSparseMatrix;
outputMat->Initialization(sparseMat->GetColNumber(),sparseMat->GetRowNumber());
nRow=sparseMat->GetRowNumber();
for(i=0;i<nRow;i++) {
GLKSparseMatrixNode* currentNode=(sparseMat->GetRowRoots())[i];
for(;currentNode!=NULL;) {
GLKSparseMatrixNode* nextNode=currentNode->pNext;
outputMat->InputElement(currentNode->nColIndex,i,currentNode->data);
currentNode=nextNode;
}
}
outputMat->InputComplete(false);
}
4a2d0fdf7This file is pinned to the selected release. Download the original to open it in your design tool.