gmres

#include <stdio.h>
#include <mpi.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>


#define m 10
#define beta 0.9

int N,flag=0,rank,nump;

void Least_sq(double *g,double **Hm)
{
    int i=0,j=0,k=0,p=0,q=0;
    double W[m+1][m+1],si=0,ci=0;
    double temp1[m+1][m],temp2[m+1];

    for (i=0;i<m;i++){
        si=(Hm[i+1][i])/sqrt(pow(Hm[i][i],2)+pow(Hm[i+1][i],2));
        ci=(Hm[i][i])/sqrt(pow(Hm[i][i],2)+pow(Hm[i+1][i],2));
        for (j=0;j<(m+1);j++){
            for (p=0; p<m; p++) {
                temp1[j][p]=0;
            }
            for (k=0;k<(m+1);k++){
                if(j==k){
                    W[j][k]=1;
                }
                else {
                    W[j][k]=0;
                }
            }
            temp2[j]=0;
        }
        W[i][i]=ci;
        W[i][i+1]=si;
        W[i+1][i]=-si;
        W[i+1][i+1]=ci;
                    for (j=0;j<(m+1);j++){
            for(k=0;k<m;k++){
                for(p=0;p<(m+1);p++){
                    temp1[j][k]+=(W[j][p]*(Hm[p][k]));
                }
            }
        }
        for (j=0; j<(m+1); j++) {
            for (k=0; k<m; k++) {
                Hm[j][k]=temp1[j][k];
            }
        }
        for (j=0;j<(m+1);j++){
            for(p=0;p<(m+1);p++){
                temp2[j]+=W[j][p]*g[p];
            }
        }
        for (j=0; j<(m+1); j++) {
            g[j]=temp2[j];
        }
    }
}

void Back_sub(double *y,double *g,double **Hm)
{
    int i=0,j=0,k=0;
    double sum=0,help[m+1];

    for (i=0;i<(m+1);i++){
        help[i]=g[i];
    }
    for (i=(m-1);i>=0;i--){
        sum=help[i];
        if (i<(m-1)){
            for (j=(i+1);j<m;j++){
                sum-=(Hm[i][j])*help[j];
            }
        }
        help[i]=sum/(Hm[i][i]);
    }
    for (i=0;i<m;i++){
        y[i]=help[i];
    }
}

void Krylov(double **A,double *uj,double **Hm, double **u_base,int rank,
        int nump,int local_N,int init,double *comm_tot)
{
    int  i=0,j=0,k=0;
    double *w,*buffer,buff,help;
    double comm_time1=0,comm_time2=0;

    w=(double*)malloc(N*sizeof(double));
    buffer=(double*)malloc(local_N*sizeof(double));
    for (j=0;j<m;j++){
                    for(k=0;k<N;k++){
                uj[k]=u_base[k][j];
        }
        for (i=0; i<local_N; i++) {
            w[i+init]=0;
            for (k=0; k<N; k++) {
                w[i+init]+=A[i][k]*uj[k];
            }
        }

        for (i=0;i<local_N; i++) {
            buffer[i]=w[i+init];
        }
        comm_time1=MPI_Wtime();
        MPI_Allgather(buffer, local_N, MPI_DOUBLE,w, local_N,
                      MPI_DOUBLE, MPI_COMM_WORLD);
        comm_time2=MPI_Wtime();
        *comm_tot += (comm_time2-comm_time1);

        for (i=0;i<=j;i++){
            for(k=0;k<N;k++){
                    uj[k]=u_base[k][i];
            }
            buff=0;
            for (k=0; k<local_N; k++) {
                buff+=w[init+k]*uj[k+init];
            }
            comm_time1=MPI_Wtime();
            MPI_Allreduce(&buff,&help,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
            comm_time2=MPI_Wtime();
            Hm[i][j]=help;
            *comm_tot+=(comm_time2-comm_time1);
                      for(k=0;k<=N;k++){
                    w[k]-=(Hm[i][j])*(uj[k]);
            }
        }

        buff=0;
        for (k=0;k<local_N;k++){
                buff+=pow(w[k+init],2);
        }
        comm_time1=MPI_Wtime();
        MPI_Allreduce(&buff,&help,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
        comm_time2=MPI_Wtime();
        Hm[j+1][j]=help;
        *comm_tot+=(comm_time2-comm_time1);
        Hm[j+1][j]=sqrt(Hm[j+1][j]);
        if (Hm[j+1][j]==0) {
            break;
        }
        if(j<(m-1)){
            for (k=0;k<N;k++){
                    u_base[k][j+1]=w[k]/(Hm[j+1][j]);
            }
        }
    }
}

void GMRES_m(double *b,double **A,double *x,int rank, int nump,int local_N,int init,int fin)
{
    int i=0,j=0,k=0,iter=0;
    double vita=0,eps=0.000000001,sum=0,buff,*buffer;
    double e[m+1],y[m],g[m+1];
    double *x0,*r0,*uj,*temp,norm=1,norm0=1;
    double **Hm,**u_base;
    double tot_time=0,init_time=0,comm_time1=0,comm_time2=0,comm_tot=0;
        init_time=MPI_Wtime();
    x0=(double*)calloc(N,sizeof(double));
    r0=(double*)malloc(N*sizeof(double));
    buffer=(double*)malloc(local_N*sizeof(double));
    uj=(double*)malloc(N*sizeof(double));
    temp=(double*)malloc(N*sizeof(double));
        Hm=(double**)(malloc((m+1)*sizeof(double*)));
    for (i=0; i<(m+1); i++) {
        Hm[i]=(double*)(malloc(m*sizeof(double)));
    }
    u_base=(double**)(malloc(N*sizeof(double*)));
    for (i=0; i<N; i++) {
        x[i]=0;
        u_base[i]=(double*)(malloc(m*sizeof(double)));
    }
    iter=1;

    while ((norm/norm0)>=eps) {
flag=iter;

        for (i=0; i<N; i++) {
            x0[i]=x[i];
            uj[i]=0;
            r0[i]=0;
            temp[i]=0;
        }
                for (i=0; i<N; i++) {
            for (j=0; j<m; j++) {
                u_base[i][j]=0;
            }
        }
        for (i=0; i<(m+1); i++) {
            e[i]=0;
            g[i]=0;
            for (j=0; j<m; j++) {
                Hm[i][j]=0;
                y[j]=0;
            }
        }
        e[0]=1;
        for (i=0; i<local_N; i++) {
            sum=0;
            for (j=0; j<N; j++) {
                sum+=A[i][j]*x0[j];
            }
            r0[i+init]=b[i+init]-sum;
        }
        vita=0;
        buff=0;
        for (i=0; i<local_N; i++) {
                vita+=pow(r0[init+i],2);
        }
        comm_time1=MPI_Wtime();
        MPI_Allreduce(&vita,&buff,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
        comm_time2=MPI_Wtime();
        comm_tot+=(comm_time2-comm_time1);
        vita=sqrt(buff);
        if (iter==1) {
            norm0=vita;
        }
        norm=vita;
        for (i=0;i<local_N; i++) {
            uj[i+init]=r0[i+init]/vita;
            buffer[i]=uj[i+init];
        }
        comm_time1=MPI_Wtime();
        MPI_Allgather(buffer, local_N, MPI_DOUBLE,uj, local_N,
                    MPI_DOUBLE, MPI_COMM_WORLD);
        comm_time2=MPI_Wtime();
        comm_tot+=(comm_time2-comm_time1);
                for (i=0; i<N; i++) {
            u_base[i][0]=uj[i];
        }
                    for (i=0; i<(m+1); i++) {
            g[i]=vita*e[i];
        }
              Krylov(A,uj,Hm,u_base,rank,nump,local_N,init,&comm_tot);
        Least_sq(g,Hm);
        Back_sub(y,g,Hm);
        for (i=0; i<local_N; i++) {
            for (k=0; k<m; k++) {
                temp[i]+=u_base[i+init][k]*y[k];
            }
        }
        for (i=0; i<local_N; i++) {
            x[i+init]=x0[i+init]+temp[i];
        }
        for (i=0;i<local_N; i++) {
            buffer[i]=x[i+init];
        }

                comm_time1=MPI_Wtime();
        MPI_Allgather(buffer, local_N, MPI_DOUBLE,x, local_N,
                      MPI_DOUBLE, MPI_COMM_WORLD);
        comm_time2=MPI_Wtime();
        comm_tot+=(comm_time2-comm_time1);
        comm_time2=comm_tot;

        iter++;
    }

    free(x0);
    free(r0);
    free(uj);
    free(temp);
    free(Hm);
    free(u_base);
    comm_time1=MPI_Wtime()-init_time;
    MPI_Allreduce(&comm_time1,&tot_time,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    MPI_Allreduce(&comm_time2,&comm_tot,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    if (rank==0) {
        printf("%d %d %.10f\n", N, nump, tot_time/nump);
           }
   }
}

void matrix_fill(double *x, double **S, double **D, double **invS)
{
    int i=0,j=0;
    double temp;

    for (i=0; i<N; i++) {
        x[i]=i+1;
        for (j=0; j<N; j++) {
            S[i][j]=0;
            D[i][j]=0;
            invS[i][j]=0;
            if (i==j) {
                S[i][j]=1;
                D[i][j]=i+1;
                invS[i][j]=1;
            }
            if ((j-i)==1) {
                S[i][j]=beta;
            }
            if (i<=j) {
                temp=(double) j-i;
                invS[i][j]=pow(-beta, temp);
            }
        }
    }

}

void matvec(double **A, double *x, double *y)
{
    double temp=0;
    int i=0,j=0;

    for (i=0; i<N; i++) {
        temp=0;
        for (j=0; j<N; j++) {
            temp+=A[i][j]*x[j];
        }
        y[i]=temp;
    }
}

void matmul(double **A, double **B, double **Y)
{
    double temp=0;
    int i=0,j=0,r=0;

    for (i=0; i<N; i++) {
        for (j=0; j<N; j++) {
            temp=0;
            for (r=0; r<N; r++) {
                temp+=A[i][r]*B[r][j];
            }
            Y[i][j]=temp;
        }
    }
}

double norm(double *x)
{
    int i=0;
    double temp,res=0;
    temp=0;
    for (i=0; i<N; i++) {
        temp+=pow(x[i], 2);
    }
    res=sqrt(temp);
    return res;
}

int main (int argc,char* argv[])
{
    MPI_Init( &argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    MPI_Comm_size(MPI_COMM_WORLD,&nump);
    double *x,*y,**A,**S,**invS,**D,**T,temp=0,help=0,eigen=0,**local_A;
    int i=0,j=0,local_N,temp_div,temp_mod,init,fin,irina;
       clock_t start1, end1 ;
    double time1,time2;
    int mem1=0,mem2=0;
      for (irina=1;irina<=16;irina++){
            N=irina*128;
    x=(double*)(malloc(N*sizeof(double)));
    y=(double*)(malloc(N*sizeof(double)));
        A=(double**)(malloc(N*sizeof(double*)));
    S=(double**)(malloc(N*sizeof(double*)));
    D=(double**)(malloc(N*sizeof(double*)));
    invS=(double**)(malloc(N*sizeof(double*)));
    T=(double**)(malloc(N*sizeof(double*)));
    for (i=0; i<N; i++) {
        S[i]=(double*)(malloc(N*sizeof(double)));
        D[i]=(double*)(malloc(N*sizeof(double)));
        invS[i]=(double*)(malloc(N*sizeof(double)));
        T[i]=(double*)(malloc(N*sizeof(double)));
    }

       matrix_fill(x,S,D,invS);
        matmul(S,D,T);
    free(D);
    free(S);
      for (i=0; i<N; i++) {
        A[i]=(double*)(malloc(N*sizeof(double)));
    }
    matmul(T,invS,A);
        matvec(A,x,y);
    free(T);
    free(invS);


       temp_mod=N%nump;
    temp_div=N/nump;
    if (temp_mod==0){
        init=rank*temp_div;
        fin=(rank+1)*temp_div-1;
    }
    else if (temp_mod!=0){
        if (rank!=(nump-1)){
            init=(rank*temp_div+rank);
            fin=((rank+1)*(temp_div+1)-1);
        }
        else if (rank==(nump-1)){
            init=(rank*temp_div+rank);
            fin=(N-1);
        }
    }

    local_N=fin-init+1;
    double *b;
    b=(double*)(malloc(N*sizeof(double)));
    for (i=0; i<N; i++) {
        x[i]=0;
        b[i]=1;
    }
    local_A=(double**)(malloc(local_N*sizeof(double*)));
    for (i=0; i<local_N; i++) {
        local_A[i]=(double*)(malloc(N*sizeof(double)));
    }
        for (i=0; i<local_N; i++) {
        for (j=0; j<N; j++) {
            local_A[i][j]=A[i+init][j];
        }
    }
    free(A);
    GMRES_m(b,local_A,x,rank,nump,local_N,init,fin);


}
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Finalize();
    return 0;
}