Matrix overflow when using PARILU

[Cielaah]

Hi,

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;
            }
        }

produces the following output:

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!

This is for the Si10H16 matrix I've downloaded from the SuiteSparse Matrix Collection.
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;
}

Any help would be greatly appreciated.

Best,
Salmah

[hartwig anzt]

Dear Salmah,

thanks for reaching out! I am happy to look into this. Would you mind sharing your problem (matrix) and maybe also the complete test routine in which you are calling the routine?

Just send it to hanzt@icl.utk.edu.

Thanks, Hartwig

[hartwig anzt]

Dear Salmah,
Sorry for the delay, but I finally found time to look into this.
It turns out that you chose a very tough problem. Even though the condition number is about 10^5, iterative methods have a hard time solving this. Also, this is not a positive definite matrix, so you won't be able to use Conjugate Gradient.
I checked with GMRES (also the matlab code) but I had no success with that solver, either. Maybe if we choose a very high restart parameter. Or a more sophisticated preconditioner than ILU(0).
Let me ask you, why are you interested in this sytem?
Thanks, Hartwig

[Cielaah]

Dear Hartwig,

thanks for looking into this. I actually want to solve a positive definite matrix, but I can't share it for data security guideline reasons. But like I said, this problem occurs with several matrices.
(Another example for this error with a positive definite matrix would be this one: https://sparse.tamu.edu/Janna/Serena)
Regardless, the error with the "overflow" should be independent of that, no?

Best,
Salmah

[hartwig anzt]

Salmah,

the ParILU preconditioner can result in overflow. This is due to the limited floating point range. I.e. dividing by a very small number (resp. 0) can result in overflow. In theory, this does not happen, but with limited precision it can. For very ill-conditioned matrices.
We are working on protecting against that, but in the meanwhile, you may have to stick to the classical ILU versions.

Sorry about that!

Hartwig

[Cielaah]

Dear Hartwig,

thank you for your response. I did not refer to an overflow for the ParILU preconditioner, but to a - in my point of view - wrong iteration over the CSR data structure. As I wrote in my original post I get the error "error: zero diagonal element in row x!". When analyzing the source of the error in more detail, I adapted the MAGMA sources by a printf and generated the output (see in my original post). The algorithm, however, processes the rows twice up to a certain point (row 4634 out of 17077 in this case), i.e. out of bounds ( this is what I initially tried to describe with "overflow"). Could you please check if this "out-of-bounds"-iteration" is intended, or if I do something wrong with the matrix set up?

Thank you, best regards,
Salmah