gmres()

void gmres

(

int

icounter

)

                        {


    int i,k,l,nr,mm,it;

    COMPLEX* r0;
    COMPLEX* b0;
    COMPLEX* v;
    COMPLEX* w;
    COMPLEX* H;
    COMPLEX* y;
    DOUBLE   eps;
    COMPLEX* cs;
    COMPLEX* sn;
    COMPLEX  sum;


    COMPLEX  gamma;
    DOUBLE   beta;
    COMPLEX* bhat;
    COMPLEX* xn;

    DOUBLE   bnrm;

    DOUBLE   error;

    int      gRows,gCols;
    int      locRows,locCols;
    int      izero = 0;
    int      ione  = 1;
    int      itemp;
    int      descY[9];
    int      descRHS[9];
    int      descsol[9];
    int      descr[9];
    int      descv[9];
    int      descw[9];
    int      info;
    int      n;
    long      m,NIT,cit;
    COMPLEX  cone  =  1.0 + I*0.0;
    COMPLEX  mcone  =  -1.0 + I*0.0;
    COMPLEX  czero  =  0.0 + I*0.0;
    DOUBLE   one  =  1.0;
    COMPLEX  tmp;
    DOUBLE   dtmp;
    double t1,t2;
    int nconv;


    DOUBLE cr,max_cr,min_cr,normr;
    int mmax,mmin,mi,d;

    char fname[255];
    FILE *f;


    gRows  = modgeominfo->ncntr + modgeominfo->nchief;
    gCols  = modgeominfo->ncntr;
    n = modgeominfo->ncntr;
    m = n;
    if(modsolinfo->restart != 0){
      m = modsolinfo->restart;
    }


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

    eps = modsolinfo->tolerance;
    H    = calloc((m+1)*(m+2),sizeof(COMPLEX));
    cs   = calloc(m+1,sizeof(COMPLEX));
    sn   = calloc(m+1,sizeof(COMPLEX));


    b0   = calloc(locCols,sizeof(COMPLEX));
    r0   = calloc(locRows,sizeof(COMPLEX));
    xn   = calloc(locRows,sizeof(COMPLEX));

    v    = calloc(locRows*(m+1),sizeof(COMPLEX));
    w    = calloc(locRows      ,sizeof(COMPLEX));


    y    = calloc(m+1,sizeof(COMPLEX));
    bhat = calloc(m+1,sizeof(COMPLEX));


    for(i =0; i<m*m;i++){
      H[i] = czero;
    }
    for(i =0; i<m;i++){
      bhat[i] = czero;
      cs[i] = czero;
      sn[i] = czero;
    }

    itemp = max(locRows,locCols);
    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); 

    
    descinit_(descsol ,&gRows ,&ione ,&gRows                  ,&ione ,&izero,&izero,&runinfo->ctxt,&gRows  ,&info);
    descinit_(descr   ,&gRows ,&ione ,&runinfo->row_block_size,&ione ,&izero,&izero,&runinfo->ctxt,&locRows,&info); 


    descinit_(descw   ,&gRows ,&ione ,&runinfo->row_block_size,&ione ,&izero,&izero,&runinfo->ctxt,&locRows,&info); 
    descinit_(descv   ,&gRows ,&ione ,&runinfo->row_block_size,&ione ,&izero,&izero,&runinfo->ctxt,&locRows,&info); 



    pdznrm2_(&gCols,&bnrm,solution->rhs,&ione,&ione,descRHS,&ione);   // beta = vec_norm2(r0,n);
    //r0 = b
    /*
    // r0 = Ax - r0
    pzgemv_("N",&gRows,&gCols,
                &cone, solution->matY  ,&ione,&ione,descY,
                       solution->sol   ,&ione,&ione,descsol,&ione,
                &cone,r0 ,&ione,&ione,descr,&ione);
    */

  /*
  printf("Ag=[");
    for(j=0;j<gRows;j++){
     for(i=0;i<gCols;i++){
       printf("%f+ %%i*%f,",creal(solution->matY[j*n+i]),cimag(solution->matY[j*n+i]));
     }
     printf(";\n");
    }
     printf("];\n");
  printf("bg=[");
    for(i=0;i<gRows;i++){
       printf("%f + %%i*%f ,\n",creal(solution->rhs[i]),cimag(solution->rhs[i]));
    }
     printf("];\n");
     */


    NIT = 1;
    if(modsolinfo->maxiterations != 0){
      NIT = modsolinfo->maxiterations;
    }

    if(rank == 0 ){
      get_filename(fname,"%s-gmres-res.out",rundetails->title);  
      f = fopen(fname,"w+");            
    }

    t1 = MPI_Wtime();
    logger(LOG_INFO,"Starting GMRES iterations: maxit=%d, restart=%d\n",NIT,m);
    cit=0;

    max_cr = .99;
    min_cr = .4;
    cr = 1.0;
    d = 3;
    mmax = modsolinfo->restart;
    mmin = 3;
    mi = m;

    // STARTING OUTER ITERATION
    for(it=0;it<NIT;it++){

      VecZeroEntries(y,(m+1));
      VecZeroEntries(v,locRows*(m+1));
      VecZeroEntries(w,locRows);
      VecZeroEntries(r0,locRows);
      VecZeroEntries(H,(m+1)*(m+2));
      // r0 = Ax - r0
      pzcopy_(&gRows,solution->rhs,&ione,&ione,descRHS,&ione,r0,&ione,&ione,descr,&ione);
      pzgemv_("N",&gRows,&gCols,
                  &mcone, solution->matY ,&ione,&ione,descY,
                          solution->sol  ,&ione,&ione,descsol,&ione,
                  &cone, r0              ,&ione,&ione,descr,&ione);

      //pzaxpy_(&gRows,&mcone,solution->rhs,&ione,&ione,descRHS,&ione,


      pdznrm2_(&gCols,&beta,r0,&ione,&ione,descr,&ione);   // beta = vec_norm2(r0,n);


      /*
         // Trying A. H. Baker, E. R. Jessup, Tz. V. Kolev "A simple strategy for varying the restart parameter in GMRES(m)"
         //   Journal of Computational and Applied Mathematics
         //    October 2007
      if(it > 0){
        dtmp = normr;
        pdznrm2_(&gCols,&normr,r0,&ione,&ione,descr,&ione);   // beta = vec_norm2(r0,n);
        cr = normr/dtmp;
      }



      if(cr > max_cr || it == 0){
        logger(LOG_INFO," (cr > max_cr, %f) changing m: %d -> %d\n",cr, m,mmax);
        m = mmax;
      }else if (cr < min_cr){
        m = m;
      }else{
        if(mi-d >=mmin){
          logger(LOG_INFO," (%f) changing m: %d -> %d\n",cr, m,m-d);
          m = m- d;
        }else{
          logger(LOG_INFO," (%f) changing m: %d -> %d\n",cr, m,mmax);
          m =mmax;
        }
      }
         */

  

  
      pzcopy_(&gRows,r0,&ione,&ione,descr,&ione, v,&ione,&ione,descv,&ione);
      tmp = one/beta;
      pzscal_(&gRows,&tmp,v,&ione,&ione,descv,&ione);// vec_mult_scal(v,r0,ONE/beta,n); v = r/beta
  
      bhat[0] = beta;
      mm = m;
  
      normr = beta;

      // GMRES ITERATION
      for(i =0; i<m;i++){
        cit++;
        if(0 == rank){
          logger(LOG_DEBUG,"Iteration #%ld/#%d\n",cit,m*NIT);
        }
  
        pzgemv_("N",&gRows,&gCols,
                    &cone, solution->matY  ,&ione,&ione,descY,
                           &(v[i*locRows]) ,&ione,&ione,descv,&ione,
                    &czero,w               ,&ione,&ione,descw,&ione); // matmult(w,A,&v[i*n],n);         // w = A*v


  
        for(k=0;k<=i;k++){
            // H[k,i] = <w,v[k]>
            // w = w - H[k,i]v[k]
            // H[k*m + i] = vec_scalar_tr(w,&v[k*n],n); // H(k,i) = w*v(k)
            // vec_sum2(w,w,&v[k*n],-H[k*m + i],n);  // w = w- H(k,i)*v(k)
  
            pzdotc_(&gRows,&tmp,w    ,&ione,&ione,descw,&ione,
                                &v[k*locRows],&ione,&ione,descv,&ione);
            H[k*m + i]  = tmp;
            tmp = -tmp;
            pzaxpy_(&gRows,&tmp,&v[k*locRows],&ione,&ione,descv,&ione,w,&ione ,&ione,descw,&ione);
  
        }
  
  
        pdznrm2_(&gRows,&dtmp,w,&ione,&ione,descw,&ione);//   H[(i+1)*m + i] = vec_norm2(w,n);          // H(i+1,i) = ||w||
        H[(i+1)*m+i] = dtmp;
  
        pzcopy_(&gRows,w,&ione,&ione,descw,&ione, &v[(i+1)*locRows],&ione,&ione,descv,&ione);
        tmp = cone/(H[(i+1)*m+i]);
        pzscal_(&gRows,&tmp,&v[(i+1)*locRows],&ione,&ione,descv,&ione); // vec_mult_scal(&v[(i+1)*n],w,ONE/(H[(i+1)*m+i]),n); // v(i+1) = w/H(i+1,i)
//        MPI_Barrier(MPI_COMM_WORLD);
  
        // Givens Rotation
        for(k=0;k<i;k++){
             gamma         =  sn[k]*H[(k+1)*m+i] + cs[k]*H[k*m+i]; // gamma = sn(k)*H(k+1,i) + cs(k)*H(k,i)
             H[(k+1)*m +i] = -sn[k]*H[(k  )*m+i] + cs[k]*H[(k+1)*m+i]; // H(k+1,i) = -sn(k)H(k,i)+cs(k)H(k+1,i)
             H[(k)*m +i]   =  gamma; // H(k,i) = gamma
        }
        rotmat(&cs[i], &sn[i], H[i*m+i], H[(i+1)*m+i]);
  

        H[i*m+i] = cs[i]*H[i*m+i] + sn[i]*H[(i+1)*m+i]; // H(i,i) = cs(i)H(i,i) + sn(i)*H(i+1,i)
        H[(i+1)*m+i] = 0.0;                           // H(i+1,i) = 0.0
  
    //for(j =0; j<locRows;j++)      { printf("%f %f\n",creal(v[locRows + j]),cimag(v[locRows + j])); }
        bhat[i+1] = -sn[i] *bhat[i];
        bhat[i]   =  cs[i] *bhat[i];
  
        error  = cabs(bhat[i+1])/bnrm;
        MPI_Bcast(&error,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
        if(rank == 0 ){
          fprintf(f,"%ld %d %d %10.8E\n",cit,it,i,error);
        }
//        MPI_Barrier(MPI_COMM_WORLD);
        if(error< eps){
          nr = i;
          nconv = cit;
          if(rank == 0){
            logger(LOG_INFO,"Reached convergence: error=%+10.8E at iteration %ld\n",error,cit);
          }
          goto SOL;
        }
  
  
      }
      --i;
      nr = m-1;
      y[nr] = bhat[nr]/H[nr*m+nr];
  
      SOL:
  
  
      for(k=i;k>=0;--k){
        sum = 0.0;
        for(l=k+1;l<=i;++l){
          sum += H[k*m+l]*y[l];
        }
        y[k] = (bhat[k] - sum) / H[k*m+k];
      }

      /*
      dtmp=0.0;
      for(l=0;l<=nr;l++){
        dtmp+=cabs(y[l])*cabs(y[l]);
      }
      printf("%f\n",sqrt(dtmp));
      */
 
      for(l=0;l<=nr;l++){
        pzaxpy_(&gRows,&y[l],&v[l*locRows],&ione, &ione, descv  , &ione,
                       solution->sol      ,&ione, &ione, descsol, &ione); //vec_sum2(x,x,&v[l*n],y[l],n); // x = x + y[l]*v
      } 


      if(error < eps){
        break;
      }
      // r0 = Ax - r0
      /*
         // Trying A. H. Baker, E. R. Jessup, Tz. V. Kolev "A simple strategy for varying the restart parameter in GMRES(m)"
         //   Journal of Computational and Applied Mathematics
         //    October 2007
      pzcopy_(&gRows,solution->rhs,&ione,&ione,descRHS,&ione,r0,&ione,&ione,descr,&ione);
      pzgemv_("N",&gRows,&gCols,
                  &mcone, solution->matY  ,&ione,&ione,descY,
                         solution->sol   ,&ione,&ione,descsol,&ione,
                  &cone, r0             ,&ione,&ione,descr,&ione);
      dtmp = normr;
      pdznrm2_(&gCols,&normr,r0,&ione,&ione,descr,&ione);   // beta = vec_norm2(r0,n);
      cr = normr/dtmp;
      */
   }



   MPI_Barrier(MPI_COMM_WORLD);
   t2 = MPI_Wtime();
   logger(LOG_DEBUG," GMRES iteration took %f s\n",(t2-t1));

   logger(LOG_INFO," GMRES: %d iterations, residual:%f\n",cit,error);
   //END CHECK
   logger(LOG_INFO,"Broadcasting solution (%d complex values)\n",gRows);
   MPI_Bcast(solution->sol,gRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
   logger(LOG_INFO," Solution sent to all nodes\n");

   if(rank == 0 ){
     fclose(f);
   }

   free(H);
   free(cs);
   free(sn);


   free(b0);
   free(r0);
   free(xn);
   free(v);
   free(w);
   free(y);
   free(bhat);

}