## Solving general banded system using ScaLAPACK from C

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

### Solving general banded system using ScaLAPACK from C

Hi there,

as the subject states, I am looking for the solution of the following instance of a banded system using ScaLAPACK:

Code: Select all
0.00,    0.00,    0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
0.00,    0.00,    0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
0.00,    0.00,   41.11,   41.11,   41.11,   41.11,   41.11,   41.11,
1.00,  -88.89,  -88.89,  -88.89,  -88.89,  -88.89,  -88.89,   -2.70,
47.78,   47.78,   47.78,   47.78,   47.78,   47.78,    4.00,    0.00,
0.00,    0.00,    0.00,    0.00,    0.00,   -1.00,    0.00,    0.00,

Notice that this matrix has already been described using LAPACK's banded data storage scheme. My concern: I have not been able to do it... nor I have been able to find a working example of a ScaLAPACK solution of a banded system from C... u.u'

These are the characteristics of the system (hopefully they are correct):

Code: Select all
/*! Set parameters: */
/* Parameters of the problem: */
the_rank = 8;
the_kl = 2;
the_ku = 1;
the_lda = the_rank;
the_ldb = the_rank;
the_nrhs = 1;
the_ja = 1;
the_ib = the_ja;

Here's my code for 1 single process:

Code: Select all
/*! Start BLACS: */
Cblacs_pinfo(&my_id, &our_np);
our_value = 0;
Cblacs_get(our_value, DEFAULT_SYSTEM_CONTEXT, &our_context);
/* Implement 1-d block-column distribution: */
our_order = 'R';
our_nprow = 1;
our_npcol = our_np;
Cblacs_gridinit(&our_context, &our_order, our_nprow, our_npcol);

/*! Set parameters: */
/* Parameters of the problem: */
the_rank = 8;
the_kl = 2;
the_ku = 1;
the_lda = the_rank;
the_ldb = the_rank;
the_nrhs = 1;
the_ja = 1;
the_ib = the_ja;

/* OPTION 2: TRANSPOSITION OF THE CORE:*/
double my_core[COEFF_PER_PROC] = {
0.00,    0.00,    0.00,    1.00,   47.78,    0.00,
0.00,    0.00,    0.00,  -88.89,   47.78,    0.00,
0.00,    0.00,   41.11,  -88.89,   47.78,    0.00,
0.00,    0.00,   41.11,  -88.89,   47.78,    0.00,
0.00,    0.00,   41.11,  -88.89,   47.78,    0.00,
0.00,    0.00,   41.11,  -88.89,   47.78,   -1.00,
0.00,    0.00,   41.11,  -88.89,    4.00,    0.00,
0.00,    0.00,   41.11,   -2.70,    0.00,    0.00
};

/* 1-d block-row distribution of the rhs: */
our_nb = the_rank/our_np;
our_mb = the_rank/our_np;
my_rhs[0] = 1.0;
my_rhs[1] = 1.0;
my_rhs[2] = 1.0;
my_rhs[3] = 1.0;
my_rhs[4] = 1.0;
my_rhs[5] = 1.0;
my_rhs[6] = 1.0;
my_rhs[7] = 0.0;

/*! Array descriptor for A: */
the_desca[0] = DTYPE_NARROW_BAND_MATRIX;  /* The descriptor type. */
the_desca[1] = our_context;               /* The BLACS context handle. */
the_desca[2] = the_rank;                  /* The numcols in global matrix. */
the_desca[3] = our_nb;                    /* The column block size. */
the_desca[4] = 0;                         /* Proc col for 1st col of mtrix. */
the_desca[5] = the_lda;                   /* Local leading dim of matrix. */
the_desca[6] = 0;                         /* Unused, reserved. */

/*! Array descriptor for B: */
the_descb[0] = DTYPE_NARROW_BAND_RHS; /* The descriptor type. */
the_descb[1] = our_context;           /* The BLACS context handle. */
the_descb[2] = the_rank;              /* The numrows in the global vector. */
the_descb[3] = our_mb;                /* The row block size. */
the_descb[4] = 0;                     /* Proc row for first row of vector. */
the_descb[5] = the_ldb;               /* Local leading dimension of vector. */
the_descb[6] = 0;                     /* Unused, reserved. */

/*! Creation of auxiliary arrays for ScaLAPACK: */
the_laf = 12*our_npcol + 3*our_nb;
the_af = (double *) malloc(the_laf*sizeof(double));
if (the_af == NULL) {
MPI_Abort(MPI_COMM_WORLD, 1);
}

the_lwork = 10*our_npcol + 4*the_nrhs;
the_work = (double *) malloc(the_lwork*sizeof(double));
if (the_work == NULL) {
MPI_Abort(MPI_COMM_WORLD, 1);
}

/*! Factorization: */
/* PDGBTRF (N, BWL, BWU, A, JA, DESCA, IPIV, AF, LAF, WORK, LWORK, INFO) */
pdgbtrf_(&the_rank, &the_kl, &the_ku, my_core, &the_ja, the_desca,
the_af, &the_laf, the_work, &the_lwork, &the_info);
if (the_info != 0) {
if (my_id == 0) {
fprintf(stderr, ":( Factorization went wrong: info = %d\n", the_info);
fprintf(stderr, "Exiting...\n");
}
MPI_Abort(MPI_COMM_WORLD, 1);
} else if (my_id == 0) {
fprintf(stdout, ":) Factorization went 0k! info = %d\n", the_info);
}

/*! Solution: */
/* PDGBTRS (TRANS, N, BWL, BWU, NRHS, A, JA, DESCA, IPIV,
B, IB, DESCB, AF, LAF, WORK, LWORK, INFO)
*/
pdgbtrs_(&the_trans, &the_rank, &the_kl, &the_ku, &the_nrhs,
my_core, &the_ja, the_desca,
my_rhs, &the_ib, the_descb, the_af, &the_laf,
the_work, &the_lwork, &the_info);
if (the_info != 0) {
if (my_id == 0) {
fprintf(stderr, ":( Factorization went wrong: info = %d\n", the_info);
fprintf(stderr, "Exiting...\n");
}
MPI_Abort(MPI_COMM_WORLD, 1);
} else if (my_id == 0) {
fprintf(stdout, ":) Solution went 0k! info = %d\n", the_info);
}

/*! Finalize BLACS: */
Cblacs_gridexit(our_context);
Cblacs_exit(EXIT_SUCCESS);

I keep getting signal 11 (seg fault) when running! What am I doing wrong? Where's the documentation or example drivers for banded systems to be solved from C?