#include "testHelper_DisconnectedMatrix.h"
#define min(x,y) ( ( (x)<(y) ) ? (x) : (y) )
#define max(x,y) ( ( (x)>(y) ) ? (x) : (y) )
/**
* Construct a random (disconnected) matrix.
*
* Input:
* symmetric : Whether the matrix should be symmetric
* dummies : Whether the matrix should contain dummy diagonal nonzeros
* colWeights : Whether the matrix should be weighted with column weights
* numCompMin : Minimum number of connected components
* numCompMax : Maximum number of connected components
*
* Output:
* pA : The generated matrix
* pNumComponents : The number of components generated
* pComponent_m : Height of each connected component
* pComponent_n : Width of each connected component
* pComponent_weights: Weight of each connected component (effective number of nonzeros or sum of column weights)
* p_i_to_I : I = p_i_to_I[cmpnnt][i] is the row index I in the large matrix which corresponds to the row index i in component cmpnnt
* p_j_to_J : J = p_j_to_J[cmpnnt][j] is the column index J in the large matrix which corresponds to the column index j in component cmpnnt
*
* Returned memory allocations (to be freed with DestructDisconnectedMatrix()):
* pA:
* pA->i O(pA->m)
* pA->j O(pA->n)
* pA->Pstart O(2+1)
* pA->dummy O(pA->m) (only if dummy == TRUE)
* pA->ColWeights O(pA->n) (only if colWeights == TRUE)
* pComponent_m: O(numComponents)
* pComponent_n: O(numComponents) (only if symmetric == FALSE)
* pComponent_weights: O(numComponents)
* p_i_to_I: O(numComponents)
* p_i_to_I[0]: O(pA->m)
* p_j_to_J: O(numComponents) (only if symmetric == FALSE)
* p_j_to_J[0]: O(pA->n) (only if symmetric == FALSE)
*/
void ConstructDisconnectedMatrix(struct sparsematrix *pA, int symmetric, int dummies, int colWeights,
long numCompMin, long numCompMax, long *pNumComponents,
long **pComponent_m, long **pComponent_n, long **pComponent_weights,
long ***p_i_to_I, long ***p_j_to_J) {
long i, j, I, J, m, n;
long cmpnnt;
/* Determine number of components (submatrices), the dimensions of each submatrix, and their numbers of nonzeros */
long numComponents = Random1(numCompMin, numCompMax);
long *component_dummies = (long *)malloc(numComponents*sizeof(long));
long *component_m = (long *)malloc(numComponents*sizeof(long));
long *component_n;
if(symmetric || dummies) {
component_n = component_m;
}
else {
component_n = (long *)malloc(numComponents*sizeof(long));
}
long *component_nnz = (long *)malloc(numComponents*sizeof(long));
if (component_dummies == NULL || component_m == NULL || component_n == NULL || component_nnz == NULL) {
printf("Error\n");
exit(1);
}
long M=0, N=0, NNZ=0, totalDummies = 0;
for(cmpnnt=0; cmpnnt<numComponents; ++cmpnnt) {
/* component_nnz[] should uniquely identify a component; i.e. all components should have different nnz.
* This helps when identifying the components in test functions.
*/
int present;
long nonDummyDiagNonZeros[2];
long full_nnz[2];
do {
/* The dimensions should at least be above 10, such that generating the
* matrices won't take too long (or, when below 5, may become impossible).
* Additionally, we should have enough freedom to ensure each component has
* a unique number of nonzeros
*/
component_m[cmpnnt] = Random1(14, 20) + Random1(numComponents, numCompMax);
if(!symmetric && !dummies) {
component_n[cmpnnt] = Random1(14, 20) + Random1(numComponents, numCompMax);
}
long maxdim = (component_m[cmpnnt]>component_n[cmpnnt]) ? component_m[cmpnnt] : component_n[cmpnnt];
/* Determine number of nonzeros */
if(symmetric) {
component_nnz[cmpnnt] = 2*maxdim-1 + Random1(maxdim/2, maxdim);
}
else {
component_nnz[cmpnnt] = 2*maxdim-1 + Random1(maxdim, 2*maxdim);
}
component_dummies[cmpnnt] = dummies ? Random1(maxdim/4, 3*maxdim/4) : 0;
component_nnz[cmpnnt] -= component_dummies[cmpnnt];
/* Check that this number of non-dummy nonzeros is unique */
present = 0;
for(i=0; i<cmpnnt; ++i) {
if(symmetric) {
nonDummyDiagNonZeros[0] = component_m[cmpnnt] - component_dummies[cmpnnt];
full_nnz[0] = (component_nnz[cmpnnt]-nonDummyDiagNonZeros[0])*2 + nonDummyDiagNonZeros[0];
nonDummyDiagNonZeros[1] = component_m[i] - component_dummies[i];
full_nnz[1] = (component_nnz[i]-nonDummyDiagNonZeros[1])*2 + nonDummyDiagNonZeros[1];
if(full_nnz[0] == full_nnz[1]) {
present = 1;
break;
}
}
else {
if(component_nnz[i] == component_nnz[cmpnnt]) {
present = 1;
break;
}
}
}
}
while(present);
/* Update totals */
M += component_m[cmpnnt];
N += component_n[cmpnnt];
NNZ += component_nnz[cmpnnt];
totalDummies += component_dummies[cmpnnt];
}
/* Allocate space for a mapping from submatrix indices (i,j) to global indices (I,J) */
long *i_to_I_alloc, *j_to_J_alloc, **i_to_I, **j_to_J;
i_to_I_alloc = (long *)malloc(M*sizeof(long));
i_to_I = (long **)malloc(numComponents*sizeof(long *));
if(symmetric || dummies) {
j_to_J_alloc = i_to_I_alloc;
j_to_J = i_to_I;
}
else {
j_to_J_alloc = (long *)malloc(N*sizeof(long));
j_to_J = (long **)malloc(numComponents*sizeof(long *));
}
if (i_to_I_alloc == NULL || i_to_I == NULL || j_to_J_alloc == NULL || j_to_J == NULL) {
printf("Error\n");
exit(1);
}
i_to_I[0] = i_to_I_alloc;
if(!symmetric && !dummies)
j_to_J[0] = j_to_J_alloc;
for(cmpnnt=1; cmpnnt<numComponents; ++cmpnnt) {
i_to_I[cmpnnt] = &i_to_I[cmpnnt-1][component_m[cmpnnt-1]];
if(!symmetric && !dummies)
j_to_J[cmpnnt] = &j_to_J[cmpnnt-1][component_n[cmpnnt-1]];
}
/* Create a permutation of the identity */
for(I=0; I<M; ++I) {
i_to_I_alloc[I] = I;
}
if(!symmetric && !dummies) {
for(J=0; J<N; ++J) {
j_to_J_alloc[J] = J;
}
}
RandomPermute(i_to_I_alloc, NULL, NULL, NULL, 0, M-1);
if(!symmetric && !dummies)
RandomPermute(j_to_J_alloc, NULL, NULL, NULL, 0, N-1);
/* Set up matrix struct */
MMSparseMatrixInit(pA);
pA->m = M;
pA->n = N;
pA->NrNzElts = NNZ + totalDummies;
pA->NrProcs = 2; /* maximum number of parts */
pA->i = (long *)malloc(pA->NrNzElts*sizeof(long));
pA->j = (long *)malloc(pA->NrNzElts*sizeof(long));
pA->Pstart = (long *)malloc((pA->NrProcs+1)*sizeof(long));
pA->RowLambda = (int *)malloc(pA->m*sizeof(int));
pA->ColLambda = (int *)malloc(pA->n*sizeof(int));
if (pA->i == NULL || pA->j == NULL || pA->Pstart == NULL || pA->RowLambda == NULL || pA->ColLambda == NULL) {
printf("Error\n");
exit(1);
}
if(colWeights)
pA->MMTypeCode[0]='W'; /* weighted matrix */
else
pA->MMTypeCode[0]='M'; /* normal matrix */
pA->MMTypeCode[1]='C'; /* coordinate scheme */
pA->MMTypeCode[2]='P'; /* pattern only */
if(symmetric)
pA->MMTypeCode[3]='S'; /* symmetric */
else
pA->MMTypeCode[3]='G'; /* general, no symmetry */
pA->NrDummies = 0;
pA->dummy = NULL;
pA->Pstart[0] = 0;
pA->Pstart[1] = pA->NrNzElts;
pA->Pstart[2] = pA->NrNzElts;
pA->NrDummies = totalDummies;
if(dummies) {
pA->dummy = (int *)malloc(pA->m*sizeof(int));
if (pA->dummy == NULL) {
printf("Error\n");
exit(1);
}
}
/* Generate each component */
long maxdim, t = 0, t2, T, nnzFilled, dummiesUsed;
for(cmpnnt=0; cmpnnt<numComponents; ++cmpnnt) {
m = component_m[cmpnnt];
n = component_n[cmpnnt];
maxdim = (m>n) ? m : n;
nnzFilled = 0;
dummiesUsed = 0;
T = t;
/* Make sure we create connected components */
if(Random1(0,1) == 0) {
/* Here, we generate a submatrix with nonzero diagonal
* and one nonzero subdiagonal. This way, each row is
* connected with the next, and hence all are connected.
*/
for(i=0; i<maxdim; ++i) {
I = i_to_I[cmpnnt][i%m];
J = j_to_J[cmpnnt][i%n];
pA->i[t] = symmetric ? max(I,J) : I;
pA->j[t] = symmetric ? min(I,J) : J;
++t;
++nnzFilled;
if(dummies) {
pA->dummy[I] = (dummiesUsed < component_dummies[cmpnnt]) ? 1 : 0;
if(pA->dummy[I] == 1) {
++dummiesUsed;
--nnzFilled;
}
}
if(i < maxdim-1) {
I = i_to_I[cmpnnt][(i+1)%m];
J = j_to_J[cmpnnt][i%n];
pA->i[t] = symmetric ? max(I,J) : I;
pA->j[t] = symmetric ? min(I,J) : J;
++t;
++nnzFilled;
}
}
}
else {
/* Here, we generate a submatrix with nonzero diagonal
* and a dense first column. This way, all rows are
* connected with the first column, hence all rows and
* columns are connected.
*/
for(i=0; i<maxdim; ++i) {
I = i_to_I[cmpnnt][i%m];
J = j_to_J[cmpnnt][i%n];
pA->i[t] = symmetric ? max(I,J) : I;
pA->j[t] = symmetric ? min(I,J) : J;
++t;
++nnzFilled;
if(dummies) {
pA->dummy[I] = (dummiesUsed < component_dummies[cmpnnt]) ? 1 : 0;
if(pA->dummy[I] == 1) {
++dummiesUsed;
--nnzFilled;
}
}
if(i > 0) {
I = i_to_I[cmpnnt][i%m];
J = j_to_J[cmpnnt][0];
pA->i[t] = symmetric ? max(I,J) : I;
pA->j[t] = symmetric ? min(I,J) : J;
++t;
++nnzFilled;
}
}
}
/* Fill the submatrix until the desired number of nonzeros */
while(nnzFilled < component_nnz[cmpnnt]) {
i = Random1(0, m-1);
j = Random1(0, n-1);
I = i_to_I[cmpnnt][i%m];
J = j_to_J[cmpnnt][j%n];
I = symmetric ? max(I,J) : I;
J = symmetric ? min(I,J) : J;
/* Make sure this element is not already nonzero.
* This search is expensive in terms of complexity, but
* as we are just testing, this does not matter very much.
*/
int present = 0;
for(t2=T; t2<t; ++t2) {
if(pA->i[t2] == I && pA->j[t2] == J) {
present = 1;
break;
}
}
if(present) {
continue;
}
pA->i[t] = I;
pA->j[t] = J;
++t;
++nnzFilled;
}
}
if(symmetric) {
/* Update nonzero count, as we now only have counted the elements in the lower-left triangle */
NNZ = 0;
long nonDummyDiagNonZeros;
for(cmpnnt=0; cmpnnt<numComponents; ++cmpnnt) {
nonDummyDiagNonZeros = component_m[cmpnnt] - component_dummies[cmpnnt];
component_nnz[cmpnnt] = (component_nnz[cmpnnt]-nonDummyDiagNonZeros)*2 + nonDummyDiagNonZeros;
NNZ += component_nnz[cmpnnt];
}
}
/* Return value pointers */
*pNumComponents = numComponents;
*pComponent_m = component_m;
*pComponent_n = component_n;
*pComponent_weights = component_nnz;
*p_i_to_I = i_to_I;
*p_j_to_J = j_to_J;
if(colWeights) {
/* Reassign weights based on columns */
long *component_weights = *pComponent_weights;
pA->NrColWeights = N;
pA->ColWeights = malloc( N * sizeof( long ) );
if (pA->ColWeights == NULL) {
printf("Error\n");
exit(1);
}
for(J=0; J<N; ++J) {
pA->ColWeights[J] = Random1(20, 50);
}
for(cmpnnt=0; cmpnnt<numComponents; ++cmpnnt) {
component_weights[cmpnnt] = 0;
for(j=0; j<component_n[cmpnnt]; ++j) {
component_weights[cmpnnt] += pA->ColWeights[j_to_J[cmpnnt][j]];
}
/* Make sure each component has a unique weight, to be able to use it as identification */
while(TRUE) {
int found = FALSE;
for(i=0; i<cmpnnt; ++i) {
if(component_weights[cmpnnt] == component_weights[i]) {
found = TRUE;
++pA->ColWeights[j_to_J[cmpnnt][0]];
++component_weights[cmpnnt];
break;
}
}
if(!found) {
break;
}
}
}
}
free(component_dummies);
} /* end ConstructDisconnectedMatrix */
/**
* Free memory for a matrix constructed with ConstructDisconnectedMatrix().
* See ConstructDisconnectedMatrix() for details on the parameters.
*/
void DestructDisconnectedMatrix(struct sparsematrix *pA, int symmetric, int dummies, int colWeights,
long **pComponent_m, long **pComponent_n, long **pComponent_weights,
long ***p_i_to_I, long ***p_j_to_J) {
free(*pComponent_weights);
free(*pComponent_m);
free(*p_i_to_I[0]);
free(*p_i_to_I);
if(!symmetric && !dummies) {
free(*pComponent_n);
free(*p_j_to_J[0]);
free(*p_j_to_J);
}
MMSparseMatrixFreeMemory(pA);
} /* end DestructDisconnectedMatrix */