#include "stdafx.h"
/*
* -- SuperLU MT routine (version 2.0) --
* Lawrence Berkeley National Lab, Univ. of California Berkeley,
* and Xerox Palo Alto Research Center.
* September 10, 2007
*
*/
#include <stdlib.h>
#include "hnum_pdsp_defs.h"
#include "hnum_cblas.h"
namespace harlinn
{
namespace numerics
{
namespace SuperLU
{
namespace Double
{
/*
* Generate a banded square matrix A, with dimension n and semi-bandwidth b.
*/
void dband(int n, int b, int nonz, double **nzval, int **rowind, int **colptr)
{
int iseed[] = {1992,1993,1994,1995};
register int i, j, ub, lb, ilow, ihigh, lasta = 0;
double *a;
int *asub, *xa;
double *val;
int *row;
printf("A banded matrix.");
dallocateA(n, nonz, nzval, rowind, colptr); /* Allocate storage */
a = *nzval;
asub = *rowind;
xa = *colptr;
ub = lb = b;
#if 0 // EH: Given the content of the iseed array, this code is just a waste of cpu cycles
for (i = 0; i < 4; ++i) iseed[i] = abs( iseed[i] ) % 4096;
if ( iseed[3] % 2 != 1 ) ++iseed[3];
#endif
for (j = 0; j < n; ++j)
{
xa[j] = lasta;
val = &a[lasta];
row = &asub[lasta];
ilow = SUPERLU_MAX(0, j - ub);
ihigh = SUPERLU_MIN(n-1, j + lb);
for (i = ilow; i <= ihigh; ++i)
{
val[i-ilow] = dlaran_(iseed);
row[i-ilow] = i;
}
lasta += ihigh - ilow + 1;
} /* for j ... */
xa[n] = lasta;
}
/*
* Generate a block diagonal matrix A.
*/
void
dblockdiag(int nb, /* number of blocks */
int bs, /* block size */
int nonz, double **nzval, int **rowind, int **colptr)
{
int iseed[] = {1992,1993,1994,1995};
register int i, j, b, n, lasta = 0, cstart, rstart;
double *a;
int *asub, *xa;
double *val;
int *row;
n = bs * nb;
printf("A block diagonal matrix: nb %d, bs %d, n %d\n", nb, bs, n);
dallocateA(n, nonz, nzval, rowind, colptr); /* Allocate storage */
a = *nzval;
asub = *rowind;
xa = *colptr;
#if 0 // EH: Given the content of the iseed array, this code is just a waste of cpu cycles
for (i = 0; i < 4; ++i) iseed[i] = abs( iseed[i] ) % 4096;
if ( iseed[3] % 2 != 1 ) ++iseed[3];
#endif
for (b = 0; b < nb; ++b)
{
cstart = b * bs; /* start of the col # of the current block */
rstart = b * bs; /* start of the row # of the current block */
for (j = cstart; j < cstart + bs; ++j)
{
xa[j] = lasta;
val = &a[lasta];
row = &asub[lasta];
for (i = 0; i < bs; ++i)
{
val[i] = dlaran_(iseed);
row[i] = i + rstart;
}
lasta += bs;
} /* for j ... */
} /* for b ... */
xa[n] = lasta;
}
};
};
};
};