Skip to content

Instantly share code, notes, and snippets.

@MrPanch

MrPanch/thomas_p

Created May 21, 2020 17:56
Show Gist options
  • Save MrPanch/d1e48544ba666c7ac3c361fe33b396e4 to your computer and use it in GitHub Desktop.
Save MrPanch/d1e48544ba666c7ac3c361fe33b396e4 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. MrPanch created this gist May 21, 2020.
    121 changes: 121 additions & 0 deletions thomas_p
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,121 @@
    void ThomasAlgorithm_P(int mynode, int numnodes,
    int N, double *b, double *a, double *c, double *x, double *q){
    int i;
    int rows_local,local_offset;
    double S[2][2],T[2][2],s1tmp,s2tmp;
    double *l,*d,*y;
    MPI_Status status;

    l = new double[N];
    d = new double[N];
    y = new double[N];

    for(i=0;i<N;i++)
    l[i] = d[i] = y[i] = 0.0;

    S[0][0] = S[1][1] = 1.0;
    S[1][0] = S[0][1] = 0.0;

    rows_local = (int) floor((double)N/numnodes);
    local_offset = mynode*rows_local;

    if(mynode==0){
    s1tmp = a[local_offset]*S[0][0];
    S[1][0] = S[0][0];
    S[1][1] = S[0][1];
    S[0][1] = a[local_offset]*S[0][1];
    S[0][0] = s1tmp;
    for(i=1;i<rows_local;i++){
    s1tmp = a[i+local_offset]*S[0][0] - b[i+local_offset-1]*c[i+local_offset-1]*S[1][0];
    s2tmp = a[i+local_offset]*S[0][1] - b[i+local_offset-1]*c[i+local_offset-1]*S[1][1];
    S[1][0] = S[0][0];
    S[1][1] = S[0][1];
    S[0][0] = s1tmp;
    S[0][1] = s2tmp;
    }
    }
    else{
    for(i=0;i<rows_local;i++){
    s1tmp = a[i+local_offset]*S[0][0] - b[i+local_offset-1]*c[i+local_offset-1]*S[1][0];
    s2tmp = a[i+local_offset]*S[0][1] - b[i+local_offset-1]*c[i+local_offset-1]*S[1][1];
    S[1][0] = S[0][0];
    S[1][1] = S[0][1];
    S[0][0] = s1tmp;
    S[0][1] = s2tmp;
    }
    }

    for(i=0; i<=log2(numnodes);i++){
    if(mynode+pow(2.0,i) < numnodes)
    MPI_Send(S,4,MPI_DOUBLE,int(mynode+pow(2.0,i)),0,MPI_COMM_WORLD);
    if(mynode-pow(2.0,i)>=0){
    MPI_Recv(T,4,MPI_DOUBLE,int(mynode-pow(2.0,i)),0,MPI_COMM_WORLD,&status);
    s1tmp = S[0][0]*T[0][0] + S[0][1]*T[1][0];
    S[0][1] = S[0][0]*T[0][1] + S[0][1]*T[1][1];
    S[0][0] = s1tmp;
    s1tmp = S[1][0]*T[0][0] + S[1][1]*T[1][0];
    S[1][1] = S[1][0]*T[0][1] + S[1][1]*T[1][1];
    S[1][0] = s1tmp;
    }
    }

    d[local_offset+rows_local-1] = (S[0][0] + S[0][1])/(S[1][0] + S[1][1]);
    if(mynode == 0){
    MPI_Send(&d[local_offset+rows_local-1],1,MPI_DOUBLE,1,0,MPI_COMM_WORLD);
    }
    else{
    MPI_Recv(&d[local_offset-1],1,MPI_DOUBLE,mynode-1,0,MPI_COMM_WORLD,&status);
    if(mynode != numnodes-1)
    MPI_Send(&d[local_offset+rows_local-1],1,MPI_DOUBLE,mynode+1,0,MPI_COMM_WORLD);
    }


    if(mynode == 0){
    l[0] = 0;
    d[0] = a[0];
    for(i=1;i<rows_local-1;i++){
    l[local_offset+i] = b[local_offset+i-1]/d[local_offset+i-1];
    d[local_offset+i] = a[local_offset+i] - l[local_offset+i]*c[local_offset+i-1];
    }
    l[local_offset+rows_local-1] = b[local_offset+rows_local-2]/d[local_offset+rows_local-2];
    }
    else{
    for(i=0;i<rows_local-1;i++){
    l[local_offset+i] = b[local_offset+i-1]/d[local_offset+i-1];
    d[local_offset+i] = a[local_offset+i] - l[local_offset+i]*c[local_offset+i-1];
    }
    l[local_offset+rows_local-1] = b[local_offset+rows_local-2]/d[local_offset+rows_local-2];
    }

    /***********************************************************************************/

    if(mynode>0)
    d[local_offset-1] = 0;

    double * tmp = new double[N];
    for(i=0;i<N;i++)
    tmp[i] = d[i];
    MPI_Allreduce(tmp,d,N,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    for(i=0;i<N;i++)
    tmp[i] = l[i];
    MPI_Allreduce(tmp,l,N,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    delete[] tmp;

    if(mynode ==0){
    /* Forward Substitution [L][y] = [q] */
    y[0] = q[0];
    for(i=1;i<N;i++)
    y[i] = q[i] - l[i]*y[i-1];

    /* Backward Substitution [U][x] = [y] */
    x[N-1] = y[N-1]/d[N-1];
    for(i=N-2;i>=0;i--)
    x[i] = (y[i] - c[i]*x[i+1])/d[i];

    }

    delete[] l;
    delete[] y;
    delete[] d;
    return;
    }