/* testR2D.c */
#include "../InpMtx.h"
#include "../../EGraph.h"
#include "../../Coords.h"
#include "../../Drand.h"
#include "../../timings.h"
/*--------------------------------------------------------------------*/
void
generateMatrix (
int eadj[],
Coords *coords,
double mtxent[],
int msglvl,
FILE *msgFile
) ;
void
loadMatrices (
EGraph *egraph,
Coords *coords,
InpMtx *inpmtx,
int msglvl,
FILE *msgFile
) ;
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------------------
test the InpMtx object
(1) read in an EGraph file that contains the vertex lists
for triangular element
(2) read in a Coords file that contains the vertex coordinates
(3) for each element
(4) create an elemental matrix and insert into InpMtx
(5) assemble the matrix
created -- 96jul05, cca
--------------------------------------------------------------
*/
int
main ( int argc, char *argv[] )
{
char *CoordsFileName, *EGraphFileName, *outFileName ;
Coords *coords ;
DenseMtx *X, *Y, *Y2 ;
double alpha[2] = {1.0, 0.0} ;
double *x, *y ;
double *mtxent ;
double t1, t2 ;
Drand *drand ;
EGraph *egraph ;
FILE *msgFile ;
int coordType, esize, ielem, msglvl, nelem, nvtx, rc, seed ;
int *colOldToNew, *eadj, *rowOldToNew ;
InpMtx *inpmtx ;
IV *colOldToNewIV, *rowOldToNewIV ;
/*
---------------
check the input
---------------
*/
if ( argc != 8 ) {
fprintf(stdout,
"\n\n usage : %s msglvl msgFile EGraphFile CoordsFile "
"\n coordType seed outFile"
"\n msglvl -- message level"
"\n msgFile -- message file"
"\n EGraphFile -- file that contains the element graph"
"\n must be *.egraphf or *.egraphb"
"\n CoordsFile -- file that contains the coordinates object"
"\n must be *.coordsf or *.coordsb"
"\n coordType -- coordinate type"
"\n 1 --> by rows"
"\n 2 --> by columns"
"\n 3 --> by chevrons"
"\n seed -- random number seed"
"\n outFile -- file to contain the InpMtx object"
"\n must be *.inpmtxb or *.inpmtxf"
"\n", argv[0]) ;
return(0) ;
}
msglvl = atoi(argv[1]) ;
if ( strcmp(argv[2], "stdout") == 0 ) {
msgFile = stdout ;
} else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
fprintf(stderr, "\n fatal error in %s"
"\n unable to open file %s\n",
argv[0], argv[2]) ;
return(-1) ;
}
EGraphFileName = argv[3] ;
CoordsFileName = argv[4] ;
coordType = atoi(argv[5]) ;
seed = atoi(argv[6]) ;
outFileName = argv[7] ;
/*
--------------
echo the input
--------------
*/
fprintf(msgFile, "\n input to %s"
"\n msglvl = %d"
"\n msgFile = %s"
"\n EGraphFile = %s"
"\n CoordsFile = %s"
"\n coordType = %d"
"\n seed = %d"
"\n",
argv[0], msglvl, argv[2], EGraphFileName, CoordsFileName,
coordType, seed) ;
fflush(msgFile) ;
/*
-------------------------
read in the EGraph object
-------------------------
*/
egraph = EGraph_new() ;
MARKTIME(t1) ;
rc = EGraph_readFromFile(egraph, EGraphFileName) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %9.5f : read in egraph from file %s",
t2 - t1, EGraphFileName) ;
if ( rc != 1 ) {
fprintf(msgFile,
"\n return value %d from EGraph_readFromFile(%p,%s)",
rc, egraph, EGraphFileName) ;
exit(-1) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n after reading EGraph object from file %s",
EGraphFileName) ;
EGraph_writeForHumanEye(egraph, msgFile) ;
fflush(msgFile) ;
}
nvtx = egraph->nvtx ;
nelem = egraph->nelem ;
/*
-------------------------
read in the Coords object
-------------------------
*/
coords = Coords_new() ;
MARKTIME(t1) ;
rc = Coords_readFromFile(coords, CoordsFileName) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %9.5f : read in egraph from file %s",
t2 - t1, CoordsFileName) ;
if ( rc != 1 ) {
fprintf(msgFile,
"\n return value %d from Coords_readFromFile(%p,%s)",
rc, coords, CoordsFileName) ;
exit(-1) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n after reading Coords object from file %s",
CoordsFileName) ;
Coords_writeForHumanEye(coords, msgFile) ;
fflush(msgFile) ;
}
/*
-------------------------------
create the random number object
-------------------------------
*/
drand = Drand_new() ;
Drand_init(drand) ;
Drand_setUniform(drand, -1.0, 1.0) ;
Drand_setSeed(drand, seed) ;
/*
-------------------------
create the InpMtx object
-------------------------
*/
MARKTIME(t1) ;
inpmtx = InpMtx_new() ;
InpMtx_init(inpmtx, coordType, INPMTX_REAL_ENTRIES,
egraph->adjIVL->tsize, 0) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %9.5f : initialize the InpMtx object",
t2 - t1) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n InpMtx after initialization") ;
InpMtx_writeForHumanEye(inpmtx, msgFile) ;
}
/*
--------------------------
fill X with random numbers
--------------------------
*/
nvtx = egraph->nvtx ;
X = DenseMtx_new() ;
DenseMtx_init(X, SPOOLES_REAL, 0, 0, nvtx, 1, 1, nvtx) ;
DenseMtx_fillRandomEntries(X, drand) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n solution vector") ;
DenseMtx_writeForHumanEye(X, msgFile) ;
}
x = DenseMtx_entries(X) ;
/*
-----------------
fill Y with zeros
-----------------
*/
Y = DenseMtx_new() ;
DenseMtx_init(Y, SPOOLES_REAL, 0, 0, nvtx, 1, 1, nvtx) ;
DenseMtx_zero(Y) ;
y = DenseMtx_entries(Y) ;
/*
---------------------------------
load the matrices and compute y[]
---------------------------------
*/
MARKTIME(t1) ;
mtxent = DVinit(9, 0.0) ;
for ( ielem = 0 ; ielem < nelem ; ielem++ ) {
IVL_listAndSize(egraph->adjIVL, ielem, &esize, &eadj) ;
if ( msglvl > 3 && msgFile != NULL ) {
fprintf(msgFile, "\n\n elemental matrix %d", ielem) ;
}
generateMatrix(eadj, coords, mtxent, msglvl, msgFile) ;
InpMtx_inputRealMatrix(inpmtx, 3, 3, 1, 3, eadj, eadj, mtxent) ;
if ( msglvl > 3 && msgFile != NULL ) {
InpMtx_writeForHumanEye(inpmtx, stdout) ;
}
y[eadj[0]] += mtxent[0] * x[eadj[0]]
+ mtxent[3] * x[eadj[1]] + mtxent[6] * x[eadj[2]] ;
y[eadj[1]] += mtxent[1] * x[eadj[0]]
+ mtxent[4] * x[eadj[1]] + mtxent[7] * x[eadj[2]] ;
y[eadj[2]] += mtxent[2] * x[eadj[0]]
+ mtxent[5] * x[eadj[1]] + mtxent[8] * x[eadj[2]] ;
}
/*
-----------------------------
change to packed storage mode
-----------------------------
*/
InpMtx_changeStorageMode(inpmtx, INPMTX_SORTED) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n InpMtx after matrices loaded") ;
InpMtx_writeForHumanEye(inpmtx, msgFile) ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n right hand side vector via elemental mvm") ;
DVfprintf(msgFile, nvtx, y) ;
}
/*
---------------------------------------------------
the matrix as constructed is singular,
uncomment this code to make node 0 a dirichlet node
and so make the matrix nonsingular
---------------------------------------------------
*/
/*
{
double *dvec ;
int ii, nent ;
int *ivec1, *ivec2 ;
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
dvec = InpMtx_dvec(inpmtx) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] == 0 ) {
if ( ivec2[ii] == 0 ) {
dvec[ii] = 1.0 ;
} else {
dvec[ii] = 0.0 ;
}
}
}
}
*/
/*
-------------------------------------------------
optionally write out the InpMtx object to a file
-------------------------------------------------
*/
if ( strcmp(outFileName, "none") != 0 ) {
InpMtx_writeToFile(inpmtx, outFileName) ;
}
/*
---------------------
change to vector mode
---------------------
*/
InpMtx_changeStorageMode(inpmtx, INPMTX_BY_VECTORS) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n InpMtx in vector mode") ;
InpMtx_writeForHumanEye(inpmtx, msgFile) ;
}
/*
---------------------
compute y[] = A * x[]
---------------------
*/
Y2 = DenseMtx_new() ;
DenseMtx_init(Y2, SPOOLES_REAL, 0, 0, nvtx, 1, 1, nvtx) ;
DenseMtx_zero(Y2) ;
MARKTIME(t1) ;
InpMtx_nonsym_mmm(inpmtx, Y2, alpha, X) ;
MARKTIME(t2) ;
fprintf(msgFile,
"\n CPU %8.3f : compute matrix-vector multiply, %.3f mflops",
t2 - t1, 2*inpmtx->nent*1.e-6/(t2 - t1)) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n right hand side vector via assembled mvm") ;
DenseMtx_writeForHumanEye(Y2, msgFile) ;
}
/*
----------------
compare Y and Y2
----------------
*/
DenseMtx_sub(Y2, Y) ;
fprintf(msgFile, "\n ||error||_max = %12.4e", DenseMtx_maxabs(Y2)) ;
/*
--------------------------------------
get random row and column permutations
--------------------------------------
*/
rowOldToNewIV = IV_new() ;
IV_init(rowOldToNewIV, nvtx, NULL) ;
IV_ramp(rowOldToNewIV, 0, 1) ;
IV_shuffle(rowOldToNewIV, seed + 1) ;
rowOldToNew = IV_entries(rowOldToNewIV) ;
colOldToNewIV = IV_new() ;
IV_init(colOldToNewIV, nvtx, NULL) ;
IV_ramp(colOldToNewIV, 0, 1) ;
IV_shuffle(colOldToNewIV, seed + 2) ;
colOldToNew = IV_entries(colOldToNewIV) ;
/*
----------------------------------
permute the InpMtx object, X and Y
----------------------------------
*/
MARKTIME(t1) ;
InpMtx_permute(inpmtx, rowOldToNew, colOldToNew) ;
DenseMtx_permuteRows(X, colOldToNewIV) ;
DenseMtx_permuteRows(Y, rowOldToNewIV) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : permute matrix", t2 - t1) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n permuted matrix") ;
InpMtx_writeForHumanEye(inpmtx, msgFile) ;
fprintf(msgFile, "\n permuted solution") ;
DenseMtx_writeForHumanEye(X, msgFile) ;
fprintf(msgFile, "\n permuted right hand side") ;
DenseMtx_writeForHumanEye(Y, msgFile) ;
}
/*
-------------------------------------------------------
compute the matrix-vector product with permuted vectors
-------------------------------------------------------
*/
DenseMtx_zero(Y2) ;
MARKTIME(t1) ;
InpMtx_nonsym_mmm(inpmtx, Y2, alpha, X) ;
MARKTIME(t2) ;
fprintf(msgFile,
"\n CPU %8.3f : compute matrix-vector multiply, %.3f mflops",
t2 - t1, 2*inpmtx->nent*1.e-6/(t2 - t1)) ;
DenseMtx_sub(Y2, Y) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n error") ;
DenseMtx_writeForHumanEye(Y2, msgFile) ;
}
fprintf(msgFile, "\n ||error||_max = %12.4e", DenseMtx_maxabs(Y2)) ;
/*
----------------
free the objects
----------------
*/
EGraph_free(egraph) ;
Coords_free(coords) ;
InpMtx_free(inpmtx) ;
Drand_free(drand) ;
DenseMtx_free(X) ;
DenseMtx_free(Y) ;
DenseMtx_free(Y2) ;
DVfree(mtxent) ;
IV_free(rowOldToNewIV) ;
IV_free(colOldToNewIV) ;
fprintf(msgFile, "\n") ;
return(1) ; }
/*--------------------------------------------------------------------*/
/*
-----------------
generate a matrix
-----------------
*/
void
generateMatrix (
int eadj[],
Coords *coords,
double mtxent[],
int msglvl,
FILE *msgFile
) {
double area, dx10, dx20, dx21, dy10, dy20, dy21 ;
double x0, x1, x2, y0, y1, y2 ;
int i, v0, v1, v2 ;
v0 = eadj[0] ;
v1 = eadj[1] ;
v2 = eadj[2] ;
x0 = Coords_value(coords, 1, v0) ;
x1 = Coords_value(coords, 1, v1) ;
x2 = Coords_value(coords, 1, v2) ;
y0 = Coords_value(coords, 2, v0) ;
y1 = Coords_value(coords, 2, v1) ;
y2 = Coords_value(coords, 2, v2) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v0, x0, y0) ;
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v1, x1, y1) ;
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v2, x2, y2) ;
fflush(msgFile) ;
}
dx10 = x1 - x0 ;
dx20 = x2 - x0 ;
dx21 = x2 - x1 ;
dy10 = y1 - y0 ;
dy20 = y2 - y0 ;
dy21 = y2 - y1 ;
area = 0.5*(dx10*dy20 - dx20*dy10) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n area = %20.12e", area) ;
fflush(msgFile) ;
}
if ( area > 1.e-12 ) {
mtxent[0] = dx21*dx21 + dy21*dy21 ;
mtxent[1] = - dx21*dx20 - dy21*dy20 ;
mtxent[2] = dx21*dx10 + dy21*dy10 ;
mtxent[3] = mtxent[1] ;
mtxent[4] = dx20*dx20 + dy20*dy20 ;
mtxent[5] = - dx20*dx10 - dy20*dy10 ;
mtxent[6] = mtxent[2] ;
mtxent[7] = mtxent[5] ;
mtxent[8] = dx10*dx10 + dy10*dy10 ;
for ( i = 0 ; i < 9 ; i++ ) {
mtxent[i] /= (4.*area) ;
}
if ( msglvl > 3 && msgFile != NULL ) {
fprintf(msgFile, "\n matrix "
"\n [ %20.12e %20.12e %20.12e ] "
"\n [ %20.12e %20.12e %20.12e ] "
"\n [ %20.12e %20.12e %20.12e ] ",
mtxent[0], mtxent[3], mtxent[6],
mtxent[1], mtxent[4], mtxent[7],
mtxent[2], mtxent[5], mtxent[8]) ;
fprintf(msgFile, "\n rowsums = [ %13.5e %13.5e %13.5e ]",
mtxent[0] + mtxent[3] + mtxent[6],
mtxent[1] + mtxent[4] + mtxent[7],
mtxent[2] + mtxent[5] + mtxent[8]) ;
fprintf(msgFile, "\n loading matrix into bag") ;
fflush(msgFile) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------
load the matrices into the bag
------------------------------
*/
void
loadMatrices (
EGraph *egraph,
Coords *coords,
InpMtx *inpmtx,
int msglvl,
FILE *msgFile
) {
double area, dx10, dx20, dx21, dy10, dy20, dy21 ;
double mtxent[9] ;
double x0, x1, x2, y0, y1, y2 ;
int esize, i, ielem, nelem, nvtx, v0, v1, v2 ;
int *eadj ;
IVL *adj ;
nvtx = egraph->nvtx ;
nelem = egraph->nelem ;
adj = egraph->adjIVL ;
/*
for ( ielem = 0 ; ielem < 10 ; ielem++ ) {
*/
for ( ielem = 0 ; ielem < nelem ; ielem++ ) {
IVL_listAndSize(adj, ielem, &esize, &eadj) ;
v0 = eadj[0] ;
v1 = eadj[1] ;
v2 = eadj[2] ;
x0 = Coords_value(coords, 1, v0) ;
x1 = Coords_value(coords, 1, v1) ;
x2 = Coords_value(coords, 1, v2) ;
y0 = Coords_value(coords, 2, v0) ;
y1 = Coords_value(coords, 2, v1) ;
y2 = Coords_value(coords, 2, v2) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v0, x0, y0) ;
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v1, x1, y1) ;
fprintf(msgFile, "\n vertex %d at (%6.3f, %6.3f)", v2, x2, y2) ;
fflush(msgFile) ;
}
dx10 = x1 - x0 ;
dx20 = x2 - x0 ;
dx21 = x2 - x1 ;
dy10 = y1 - y0 ;
dy20 = y2 - y0 ;
dy21 = y2 - y1 ;
area = 0.5*(dx10*dy20 - dx20*dy10) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n area = %20.12e", area) ;
fflush(msgFile) ;
}
if ( area > 1.e-12 ) {
mtxent[0] = dx21*dx21 + dy21*dy21 ;
mtxent[1] = - dx21*dx20 - dy21*dy20 ;
mtxent[2] = dx21*dx10 + dy21*dy10 ;
mtxent[3] = mtxent[1] ;
mtxent[4] = dx20*dx20 + dy20*dy20 ;
mtxent[5] = - dx20*dx10 - dy20*dy10 ;
mtxent[6] = mtxent[2] ;
mtxent[7] = mtxent[5] ;
mtxent[8] = dx10*dx10 + dy10*dy10 ;
for ( i = 0 ; i < 9 ; i++ ) {
mtxent[i] /= (4.*area) ;
}
if ( msglvl > 3 && msgFile != NULL ) {
fprintf(msgFile, "\n matrix %d"
"\n [ %20.12e %20.12e %20.12e ] "
"\n [ %20.12e %20.12e %20.12e ] "
"\n [ %20.12e %20.12e %20.12e ] ",
ielem,
mtxent[0], mtxent[3], mtxent[6],
mtxent[1], mtxent[4], mtxent[7],
mtxent[2], mtxent[5], mtxent[8]) ;
fprintf(msgFile, "\n rowsums = [ %13.5e %13.5e %13.5e ]",
mtxent[0] + mtxent[3] + mtxent[6],
mtxent[1] + mtxent[4] + mtxent[7],
mtxent[2] + mtxent[5] + mtxent[8]) ;
fprintf(msgFile, "\n loading matrix into bag") ;
fflush(msgFile) ;
}
InpMtx_inputRealMatrix(inpmtx, 3, 3, 1, 3, eadj, eadj, mtxent) ;
if ( msglvl > 3 && msgFile != NULL ) {
InpMtx_writeForHumanEye(inpmtx, stdout) ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1