I wrote a little application that reads a sparse matrix and an optinal rhs and analyses them with user-input solver and preconditioner parameters. All works well except when I want to use PARILU as preconditioner. When printing the solver information I get an "error: zero diagonal element in row x!" followed by a floating point exception and the program stopping. Also, the x changes each time a run the program. I have tested several matrices and located the error to be in magma_djacobisetup_diagscal (line 278ff) in djacobi.cpp.
Adding a printf in magma_djacobisetup_diagscal:
Code: Select all
for( magma_int_t rowindex=0; rowindex<B.num_rows; rowindex++ ) {
magma_int_t start = (B.drow[rowindex]);
magma_int_t end = (B.drow[rowindex+1]);
for( i=start; i<end; i++ ) {
if ( B.dcol[i]==rowindex ) {
diag.val[rowindex] = 1.0/B.val[i];
printf("%f at row %d\n", diag.val[rowindex], rowindex);
break;
}
}
Code: Select all
...
1.000000 at row 17072
1.000000 at row 17073
1.000000 at row 17074
1.000000 at row 17075
1.000000 at row 17076
0.062949 at row 0
0.062974 at row 1
0.065418 at row 2
0.065445 at row 3
0.065459 at row 4
...
0.000000 at row 4629
-0.000000 at row 4630
-0.000000 at row 4631
-nan at row 4632
0.000000 at row 4633
-0.000000 at row 4634
error: zero diagonal element in row 4634!
So it seems there is some kind of overflow somewhere, but I couldn't manage to find the exact source.
Another problem, which may be unrelated, is a "solver class not supported" error each time I try to use ILU or PARILU as preconditioner with a P-solver like PCG. Is there a specific configuration to use these preconditioners or am I maybe missing something obvious?
My program code looks like this:
Code: Select all
int main( int argc, char** argv )
{
int i, precond_levels=0, sweeps=0, maxiter=1000, verbose=100;
double *rhs, atol=1e-16;
// Initialize MAGMA and create some LA structures.
magma_init();
magma_queue_t queue;
magma_queue_create( 0, &queue );
magma_d_matrix A={Magma_CSR};
// Read matrix from input.
magma_d_csr_mtx(&A, argv[1], queue);
int m = A.num_rows;
rhs = (double*) calloc(m, sizeof(double));
//DOUBLE PRECISION EVALUATION
printf("\n===========DOUBLE PRECISION EVALUATION===========\n");
magma_dopts opts;
magma_d_matrix dA={Magma_CSR};
magma_d_matrix b={Magma_CSR}, db={Magma_CSR};
magma_d_matrix x={Magma_CSR}, dx={Magma_CSR};
// If rhs is not given by user input, set x=1.
if (argc < 3) {
// Set x = 1.
for (i = 0; i < m; ++i) {
rhs[i] = 1.0;
}
} else {
FILE * file = fopen (argv[2], "r");
if(!file)
perror("file not found\n");
char line[256];
i = 0;
while(!feof(file))
{
fscanf(file, "%s", line);
double val = atof(line);
rhs[i++] = val;
}
}
// Pass the system to MAGMA.
magma_dvset( m, 1, rhs, &b, queue );
// Read solver and preconditioner configurations from file.
FILE * config = fopen ("solver_config.txt", "r");
char key[256], value[256];
if (config == NULL)
perror("file not found\n");
while (!feof(config)) {
fscanf(config, "%s = %s", key, value);
// Set parameters from file.
if (strcmp(key, "solver") == 0)
opts.solver_par.solver = strtosolver(value);
if (strcmp(key, "precond") == 0)
opts.precond_par.solver = strtoprecond(value);
if (strcmp(key, "maxiter") == 0)
maxiter = atoi(value);
if (strcmp(key, "verbose") == 0)
verbose = atoi(value);
if (strcmp(key, "levels") == 0)
precond_levels = atoi(value);
if (strcmp(key, "sweeps") == 0)
sweeps = atoi(value);
if (strcmp(key, "atol") == 0)
atol = atof(value);
}
fclose(config);
opts.solver_par.restart = 8;
opts.solver_par.maxiter = maxiter;
opts.solver_par.rtol = 1e-10;
opts.solver_par.verbose = verbose;
opts.precond_par.levels = precond_levels;
opts.precond_par.sweeps = sweeps;
opts.precond_par.atol = atol;
// Initialize the solver.
magma_dsolverinfo_init( &opts.solver_par, &opts.precond_par, queue );
// Copy the system to the device
magma_dmtransfer( A, &dA, Magma_CPU, Magma_DEV, queue );
magma_dmtransfer( b, &db, Magma_CPU, Magma_DEV, queue );
magma_dmtransfer( x, &dx, Magma_CPU, Magma_DEV, queue );
// Initialize an initial guess for the iteration vector.
magma_dvinit( &dx, Magma_DEV, b.num_rows, b.num_cols, 0.0, queue );
// Generate the preconditioner.
magma_d_precondsetup( dA, db, &opts.solver_par, &opts.precond_par, queue );
// Solving the problem.
magma_d_solver( dA, db, &dx, &opts, queue );
// Print solving information.
magma_dsolverinfo( &opts.solver_par, &opts.precond_par, queue );
printf("\n");
// Free memory.
magma_dsolverinfo_free( &opts.solver_par, &opts.precond_par, queue);
magma_dmfree( &dx, queue );
magma_dmfree( &db, queue );
magma_dmfree( &dA, queue );
magma_dmfree( &x, queue );
magma_dmfree( &b, queue );
magma_dmfree( &A, queue );
magma_queue_destroy(queue);
magma_finalize();
return 0;
}
Best,
Salmah