matyfd()

void matyfd	(	COMPLEX	cfreq,
		int	icounter
	)

Evaluates Matrix Y and RHS of the linear System. Valid for Neumann-type & Dirichlet-type Boundary Conditions. Coefficients are read from file.

Parameters

cfreq	COMPLEX frequency
icounter	index of frequency

                                       {

  int icntr;
  int ielmb,ielmu;
  int nphiu;
  int nphib;

  int mpi_ncntr,mpi_nelmu;
  COMPLEX cexpc;
  DOUBLE bb,cc,dd,thetac,estar;
  COMPLEX matval;
  COMPLEX yy = 0.0 + I*0.0;;
  int cstride,bstride;
  struct panel4* mpi_elements;

  //DOUBLE *C;

   Vec psiuvec;
   Vec psincvec;
   Vec phincvec;
   int gCRows;
   int nn;
   int mm;
   int descpsiu[9];
   int descgpsiu[9];
   int ione = 1;
   int info;


  int gRows,gCols;     // global dim of coefficient matrix
  int locRows,locCols; // local dim of coefficient matrix
  int izero = 0;
  int i,j;
   COMPLEX* ctmp;

#ifdef _ACOUSTO_DEBUG
  FILE* f;
  FILE* cmatfile;
  char fname[15];
  int ipsiu;
#endif

   FILE* fccoef = NULL;
   FILE* fbcoef = NULL;
   int nrow,nrowi,iloop,nloop,irow,nio,nrlast;
   char fcname[MAX_PATH];
   char fbname[MAX_PATH];

   int ic,is,ib;
   int csymmoffset,bsymmoffset;
   COMPLEX lambda;
   COMPLEX gamma;
   COMPLEX func;
   COMPLEX gamma_over_lambda;
   COMPLEX lambda_over_gamma;
   DOUBLE gdotn;
   struct vector g;
   struct vector un;

   MPI_Group* mpi_row_groups;
   MPI_Group mpi_group_world;
   MPI_Comm* mpi_comm;
   int* row_ranks;
   int ierr;

  // number of control points,phib and phiu points
  mpi_elements = geometry->elements;
  mpi_ncntr  = modgeominfo->ncntr + modgeominfo->nchief;
  nphiu      = modgeominfo->ncntr;
  nphib      = modgeominfo->nelmb;
  mpi_nelmu  = nphiu;

  gRows  = mpi_ncntr;
  gCols  = nphiu;


#ifdef _ACOUSTO_DEBUG
  cmatfile = fopen("CD.matrix.out","w");
  fprintf(cmatfile,"CD=[");
#endif

  logger(LOG_INFO,"bc %d\n",modsolinfo->dirneu);


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

  // nrow partition on node to limit memory usage (see manual)
  nrow  = runinfo->krow;
  if(runinfo->krow <= 0 || runinfo->krow >= locRows) nrow = locRows;

  // nrow partition on node to limit memory usage (see manual)
  if (solution->matY == NULL) {
    solution->matY = calloc(locRows*locCols,sizeof(COMPLEX)); // allocation
  } else {
    for(i = 0; i < locRows*locCols; i++)
      solution->matY[i] = COMPLEX_ZERO;
  }

  gCRows     = nphiu;
  psiuvec  = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  psincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  phincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  mm = gCRows;// * runinfo->nprows;
  nn = runinfo->npcols;
  descinit_(descpsiu   , &gCRows  ,&ione, &runinfo->row_block_size  ,&ione,&izero,&izero,&runinfo->ctxt,&gCRows  ,&info);
  descinit_(descgpsiu  , &mm      ,&ione, &gCRows                   ,&ione,&izero,&izero,&runinfo->ctxt,&gCRows  ,&info);
  pzgemr2d_(&gCRows,&ione,solution->vecPsiu ,&ione, &ione,descpsiu, psiuvec ,&ione,&ione,descgpsiu,&runinfo->ctxt);
  pzgemr2d_(&gCRows,&ione,solution->vecPsinc,&ione, &ione,descpsiu, psincvec,&ione,&ione,descgpsiu,&runinfo->ctxt);
  pzgemr2d_(&gCRows,&ione,solution->vecPhinc,&ione, &ione,descpsiu, phincvec,&ione,&ione,descgpsiu,&runinfo->ctxt);

  MPI_Bcast(psiuvec ,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
  MPI_Bcast(psincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
  MPI_Bcast(phincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);

  CREATE_VEC(solution->rhs,locRows); // allocation
  VecZeroEntries(solution->rhs,locRows); // reset

  CREATE_VEC(ctmp,locRows); // allocation


  if(runinfo->krow > 0){
    // C coefficients (calculated in ac_body_coefs.c)
    bodycoefs->C = calloc(4*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));
    bodycoefs->B = calloc(2*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));
  }

  j12(locRows,nrow,&nrlast,&nloop); // calculating last number of rows and number of loops

  if(runinfo->krow > 0){
   get_filename(fbname,"bbcoefs.bin.%d",rank); // B coefficients file name
   fbcoef = fopen(fbname,"rb");
  }
  if(runinfo->krow > 0){
    get_filename(fcname,"bccoefs.bin.%d",rank); // name of file with C Coefficients
    fccoef = fopen(fcname,"rb");
  }

  for(iloop=0;iloop<nloop;iloop++){  // LOOP ON NLOOP
   nrowi = nrow;
   if(iloop == nloop-1) nrowi = nrlast;
   cstride = nrowi*locCols;
   bstride = nrowi*locCols; // B stride

   if(runinfo->krow > 0){
     nio = fread(bodycoefs->C,sizeof(DOUBLE),4*nrowi*locCols*runinfo->nsymm,fccoef);  // reading from file
     nio = fread(bodycoefs->B,sizeof(DOUBLE),2*nrowi*locCols*runinfo->nsymm,fbcoef);  // Reading from file nrowi rows for this partition
   }
   for(irow=0;irow<nrowi;irow++){    // LOOP ON NROW
    i= j21(irow,iloop,nrow);
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);

    if(modsolinfo->knw == 1){
      if(runinfo->mycol == 0){
          if (modsolinfo->dirneu == 1) solution->rhs[i] =  -phincvec[icntr]; // DIRICHLET-TYPE
          if (modsolinfo->dirneu == 2) solution->rhs[i] =   phincvec[icntr]; // NEUMANN-TYPE
      }else{
        solution->rhs[i] = COMPLEX_ZERO;
      }
    }else{
      solution->rhs[i] = COMPLEX_ZERO;
    }

    for(j=0;j<locCols;j++){ // LOOP ON NELMB
      sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ielmu);
      ic = irow*locCols +j;
      ib = irow*locCols +j;

      estar = (icntr == ielmu ) ? 0.5 : 0.0;
      solution->matY[locRows*j + i ] += estar;
       for(is=0;is<runinfo->nsymm;is++){

        ielmb = ielmu + is*mpi_nelmu;

        csymmoffset = is * 4*nrowi*locCols;
        bsymmoffset = is * 2*nrowi*locCols;

        bb     = bodycoefs->B[ib             + bsymmoffset];
        cc     = bodycoefs->C[ic             + csymmoffset];
        dd     = bodycoefs->C[ic + cstride   + csymmoffset];
        thetac = bodycoefs->C[ic + cstride*2 + csymmoffset];
        cexpc  = CEXP(-cfreq*thetac);
  
                vec_copy(&un,mpi_elements[ielmb].n);
                lambda = bc.lambda[ielmu];
                gamma  = bc.gamma[ielmu];
                func   = bc.func[ielmu];
                gamma_over_lambda = gamma/lambda;
                lambda_over_gamma = lambda/gamma;
                g.x    = bc.g.x;
                g.y    = bc.g.y;
                g.z    = bc.g.z;
                gdotn  = vec_dot(g,un);
                matval =  (cc+cfreq*dd) * cexpc ;
                //
                // THIS IS CLASSICAL YMAT=E*I-(C+s*D). NEEDED TO BUILD NEUMANN & DIRICHLET RHSs.
                // NOTE: matval includes delay
                //
                solution->matY[locRows*j + i ] -= matval;
                yy     = solution->matY[locRows*j + i ];

                if(modsolinfo->knw == 1){
                        if (modsolinfo->dirneu == 1) solution->rhs[i] += yy * (func+gdotn)/gamma;                       // DIRICHLET
                        if (modsolinfo->dirneu == 2) solution->rhs[i] += bb * cexpc * (func+gdotn)/lambda;      // NEUMANN
                }else if(modsolinfo->knw == 2){
                        if (modsolinfo->dirneu == 1) solution->rhs[i] += yy * ((func+gdotn)/gamma - psincvec[ielmu]*lambda_over_gamma - phincvec[ielmu]);                       // DIRICHLET
                        if (modsolinfo->dirneu == 2) solution->rhs[i] += bb * cexpc*((func+gdotn)/lambda - phincvec[ielmu]*gamma_over_lambda - psincvec[ielmu]);        // NEUMANN
                }
                //
                // HERE SYSTEM MATRIX IS BUILT ON YMAT.
                //
                if (modsolinfo->dirneu == 1) solution->matY[locRows*j + i ]  = bb*cexpc + yy * lambda_over_gamma; //DIRICHLET (matrix entry is rewritten)
                if (modsolinfo->dirneu == 2) solution->matY[locRows*j + i ] += bb * gamma_over_lambda *cexpc; //NEUMANN   (matrix entry is updated)

#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,"%f + %%i*%f ,",CREAL(solution->matY[locRows*j + i ]),CIMAG(solution->matY[locRows*j + i ]));;
#endif

      }
    } // END OF LOOP ON NELMB
#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,";\n");
#endif
   } // END OF LOOP ON NROW
  } // END OF LOOP ON NLOOP

#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,"];\n");
#endif

   
  if(runinfo->krow > 0){
    free(bodycoefs->B);
    free(bodycoefs->C);
  }

  // Reducing RHS if cols in proc grid are > 1
  if(runinfo->npcols > 1){
    logger(LOG_INFO," Reducing RHS on columns\n");
    mpi_row_groups = calloc(sizeof(MPI_Group) , runinfo->nprows);
    mpi_comm = calloc(sizeof(MPI_Comm) , runinfo->nprows);
    row_ranks = calloc(sizeof(int) , runinfo->npcols);
    MPI_Comm_group(MPI_COMM_WORLD, &mpi_group_world);
    for(i=0;i<runinfo->nprows;i++){
       for(j=0;j<runinfo->npcols;j++){
         row_ranks[j] = j + i*runinfo->npcols;
       }
       ierr = MPI_Group_incl(mpi_group_world,runinfo->npcols,row_ranks,&mpi_row_groups[i]);
       ierr = MPI_Comm_create(MPI_COMM_WORLD, mpi_row_groups[i], &mpi_comm[i]);
    }
    
    for(i=0;i<runinfo->nprows;i++){
    
      if(MPI_COMM_NULL != mpi_comm[i]){
        MPI_Allreduce((void*)solution->rhs,
                      (void*)&ctmp[0],
                      locRows,
                      mpi_complex,
                      mpi_op_complex_sum,
                      mpi_comm[i]);
      }
//                    MPI_COMM_WORLD);
    }
    
    for(i=0;i<locRows;i++){
      solution->rhs[i] = ctmp[i];
//      printf("%d %f %f %f %f\n",i, CREAL(solution->rhs[i]),CIMAG(solution->rhs[i]),CREAL(phincvec[i]),CIMAG(phincvec[i]));
    }
  }
#ifdef _ACOUSTO_DEBUG
  sprintf(fname,"maty.%d%d",runinfo->myrow,runinfo->mycol);
  f = fopen(fname,"w");
  for(i=0;i<locRows;i++){    // LOOP ON ROWS
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    for(j=0;j<locCols;j++){ // LOOP on COLS
      sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ipsiu);
      fprintf(f,"%d %d %f %f\n",icntr,ipsiu,CREAL(solution->matY[i+locRows*j]), CIMAG(solution->matY[i+locRows*j]));
    }
  }
  fclose(f);

  fclose(cmatfile);
#endif
  if(psiuvec  ) free(psiuvec);
  if(psincvec ) free(psincvec);
  if(phincvec ) free(phincvec);
  if(ctmp ) free(ctmp);

}