solve_linear_system()

void solve_linear_system

(

int

icounter

)

Solves the linear system of equation stored in solution_struct. The function uses the ScaLAPACK functions to solve the system.

Parameters

[in]

icounter

                                      {

#ifdef _ACOUSTO_DEBUG
  FILE* rhsfile;
  char fname[15];
  int ipsiu;
  int i;
#endif
  int descY[9];
  int descRHS[9];

  COMPLEX* work;
  int lwork;
  int info;

  int gRows,gCols;
  int locRows,locCols;

  int descsol[9];

  int itemp;
  int izero = 0;
  int ione =1;

  int nn,mm;
  double t1,t2;



#ifdef _ACOUSTO_DEBUG
  rhsfile = fopen("rhs.sc","w");
  fprintf(rhsfile,"b=[");
  for(i=0;i<modgeominfo->ncntr;i++){
    fprintf(rhsfile,"%f + %%i*%f ,",CREAL(solution->rhs[ i ]),CIMAG(solution->rhs[ i ]));
  }
  fprintf(rhsfile,"]");
  fclose(rhsfile);
#endif

  gRows  = modgeominfo->ncntr + modgeominfo->nchief;
  gCols  = modgeominfo->ncntr;

  locRows  = numroc_(&gRows ,&runinfo->row_block_size,&runinfo->myrow,&izero,&runinfo->nprows);
  locCols  = numroc_(&gCols ,&runinfo->col_block_size,&runinfo->mycol,&izero,&runinfo->npcols);


  itemp = max(1,locRows);
  descinit_(descY   ,&gRows,&gCols,&runinfo->row_block_size,&runinfo->col_block_size,&izero,&izero,&runinfo->ctxt,&locRows,&info); 
  descinit_(descRHS ,&gRows,&ione ,&runinfo->row_block_size,&ione                        ,&izero,&izero,&runinfo->ctxt,&locRows,&info); 

  work = calloc(2,sizeof(COMPLEX));
  lwork = -1;

  pzgels_( "N", &gRows,&gCols,&ione, solution->matY, &ione,&ione,descY, solution->rhs , &ione,&ione,descRHS, work,&lwork,&info);

  lwork = (int) work[0];
 
  free(work);
  work = calloc(lwork,sizeof(COMPLEX));

  t1 = MPI_Wtime();
  pzgels_( "N", &gRows,&gCols,&ione, solution->matY, &ione,&ione,descY, solution->rhs , &ione,&ione,descRHS, work,&lwork,&info);
  t2 = MPI_Wtime();
  logger(LOG_DEBUG," Pseudoinverse evaluation took %f s\n",(t2-t1));

  free(work);

   // Collecting solution on node zero and broadcasting to everyone
   mm = gRows;// * runinfo->nprows;
   nn = runinfo->npcols;
   descinit_(descsol , &mm   ,&ione, &gRows  ,&ione,&izero,&izero,&runinfo->ctxt,&gRows  ,&info);
   pzgemr2d_(&gRows,&ione,solution->rhs,&ione, &ione,descRHS, solution->sol,&ione,&ione,descsol,&runinfo->ctxt);


   MPI_Bcast(solution->sol,gRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);


}