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

Postby ejspeiro » Thu Feb 21, 2013 3:27 pm

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' :roll:

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!
ejspeiro
 
Posts: 12
Joined: Wed Sep 26, 2012 11:30 am
Location: San Diego, CA

Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 1 guest