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?Thanks in advanced!
