银河里的星星

落在人间

日志

关于我

星星

文章分类

线性方程组求解-高斯消元

2009-04-22 21:45:43| 分类：高性能计算 | 标签： |举报 |字号大中小订阅

下载LOFTER 我的照片书 |

4    5
1.000000    4.000000    -2.000000    3.000000    6.000000
2.000000    2.000000    0.000000    4.000000    2.000000
3.000000    0.000000    -1.000000    2.000000    1.000000
1.000000    2.000000    2.000000    -3.000000    8.000000

一.基本实现

#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;
void environment_init();
typedef double data_t;
void cblas_factor(data_t *vec,int size,data_t factor);
void cblas_sub(data_t *minuend,data_t *subtrahend,int size);
data_t * cblas_max(data_t *vec,int size);
//N行，M列，M = N+1
int N,M;
//X的系数矩阵
data_t * matrix_A;
//y向量
data_t * B;
//x的解向量
data_t * X;
//通过A(i,j)得到，x系数矩阵中[i,j]元素的地址
data_t * A(int i,int j){
    return matrix_A+i*N+j;
}
void print_matrix(){
    printf("matrix:\n");
    for(int i = 0 ; i < N ; i++){
        for(int j = 0 ; j < N ; j++){
                printf("%lf ",*A(i,j));
        }
        printf("%lf\n",*(B+i));
    }

}
int main()
{
    FILE * f = fopen("data.txt","r");
    fscanf(f,"%d%d",&N,&M);
    matrix_A = (data_t *)malloc(sizeof(data_t)*N*N);
    B = (data_t *)malloc(sizeof(data_t)*N);
    X = (data_t *)malloc(sizeof(data_t)*N);
    int *shift_table = (int*)malloc(sizeof(int)*N);

    //get matrix from file
    for(int i = 0 ; i < N ; i++){
        for(int j = 0 ; j < N ; j++){
            fscanf(f,"%lf",A(i,j));
            printf("%lf ",*A(i,j));
        }
        fscanf(f,"%lf",B+i);
        printf("%lf\n",*(B+i));
    }
    //init shift table
    for(int i = 0 ; i < N ; i++){
        shift_table[i] = i;
    }
    //computer
    for(int i = 0 ; i < N-1 ; i++){
        //选主元
        int index = cblas_max(A(i,i),N-i)-A(i,0);
        //列交换
        if(index != i){

    //错误！！！
           shift_table[i] = index;
           shift_table[index] = i;

           for(int row = 0 ; row < N ; row++){
               data_t temp = *A(row,i);
               *(A(row,i)) = *A(row,index);
               *A(row,index) = temp;
           }
        }
        //主元归一化
        // 消元

    //错误！！！ *A(row,i)!= 0会让消元提前结束
        for(int row = i+1 ; *A(row,i)!= 0 && row < N ; row++){
            data_t factor = *A(i,i)/ *A(row,i);
            cblas_factor(A(row,i),N-i,factor);
            cblas_sub(A(row,i),A(i,i),N-i);
            *(B+row) *= factor;
            *(B+row) = *(B+row)-*(B+i);
        }
        print_matrix();

    }
    //回代，i代表了Xi,col代表了列号，用来计算前面已经计算的xi值的sum
    for(int i = N-1 ; i >= 0 ; i--){
        data_t sum = 0;
        for(int col = i+1 ; col < N ; col++){
            sum += *A(i,col)*X[col];
        }
        X[i] = (*(B+i) - sum) / *A(i,i);
    }
    //恢复,其中shift table保存了原来的列号
    for(int i = 0 ; i < N ; i++){
        printf("x%d:%lf ",shift_table[i],X[i]);
    }

}

//向量操作，对所有的数*factor
void cblas_factor(data_t *vec,int size,data_t factor){
    for(int i = 0 ; i < size ; i++){
        *(vec+i) = *(vec+i)*factor;
    }
}
//两个向量的减
void cblas_sub(data_t *minuend,data_t *subtrahend,int size){
    for(int i = 0 ; i < size ; i++){
        *(minuend+i) = *(minuend+i) - *(subtrahend+i);
    }

}
//寻找向量中的最大值
data_t * cblas_max(data_t *vec,int size){
    data_t * index = vec;
    data_t max = *vec;
    for(int i = 0 ; i < size ; i++){
        if(*(vec+i) > max){
            max = *(vec+i);
            index = vec+i;
        }
    }
    return index;

}

二.使用blas上的cell实现

blas参考文档：http://www.netlib.org/blas/blast-forum/blas-report.pdf

基本blas实现：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include "malloc_align.h"

void environment_init();
typedef double data_t;
int N,M;//N行，M列，M = N+1
data_t * matrix_A;//X的系数矩阵
data_t * B;//y向量
data_t * X;//x的解向量
int *shift_table;

//通过A(i,j)得到，x系数矩阵中[i,j]元素的地址
data_t * A(int i,int j){
    return matrix_A+i*N+j;
}
void print_matrix();

int main()
{
    FILE * f = fopen("data.txt","r");
    int i = 0,j = 0,index = 0,row = 0,temp = 0;

    fscanf(f,"%d%d",&N,&M);
    matrix_A = (data_t *)_malloc_align(sizeof(data_t)*N*N,7);
    B = (data_t *)_malloc_align(sizeof(data_t)*N,7);
    X = (data_t *)_malloc_align(sizeof(data_t)*N,7);
    shift_table = (int*)_malloc_align(sizeof(int)*N,7);

    //get matrix from file
    for(i = 0 ; i < N ; i++){
        for(j = 0 ; j < N ; j++){
            fscanf(f,"%lf",A(i,j));
            printf("%lf ",*A(i,j));
        }
        fscanf(f,"%lf",B+i);
        printf("%lf\n",*(B+i));
    }
    //init shift table
    for(i = 0 ; i < N ; i++){
        shift_table[i] = i;
    }
    //computer
    for(i = 0 ; i < N-1 ; i++){
        //选主元
    index = cblas_idamax(N-i,A(i,i),1)+i;//cblas_idamax 选择最大绝对值的blas库
    printf("index:%d\n",index);
        //列交换
        if(index != i){
           temp = shift_table[i];
       shift_table[i] = shift_table[index];
           shift_table[index] = temp;
           for(row = 0 ; row < N ; row++){
               data_t temp = *A(row,i);
               *(A(row,i)) = *A(row,index);
               *A(row,index) = temp;
           }
        }
    for(row = 0 ; row < N ; row++)
        printf("%d ",shift_table[row]);
    printf("\n");
        //主元归一化
        // 消元
        for(row = i+1 ; row < N ; row++){
        if(*A(row,i) == 0) continue;
        data_t factor = -1*(*A(row,i)/ *A(i,i));
        cblas_daxpy(N-i,factor,A(i,i),1,A(row,i),1);//cblas_daxpy( N, alpha,X,incX, Y, incY); Y=Y+alpha*X
        *(B+row) = *(B+row) + *(B+i)*factor;
        }
        print_matrix();

    }
    //回代，i代表了Xi,col代表了列号，用来计算前面已经计算的xi值的sum
    for(i = N-1 ; i >= 0 ; i--){
        data_t sum = 0;
        for( j = i+1 ; j < N ; j++){
            sum += *A(i,j)*X[j];
        }
        X[i] = (*(B+i) - sum) / *A(i,i);
    }
    //恢复,其中shift table保存了原来的列号
    for(i = 0 ; i < N ; i++){
        printf("x%d:%lf ",shift_table[i],X[i]);
    }


}
void print_matrix(){
    int i,j;
    printf("matrix:\n");
    for(i = 0 ; i < N ; i++){
        for(j = 0 ; j < N ; j++){
                printf("%lf ",*A(i,j));
        }
        printf("%lf\n",*(B+i));
    }

}

三.使用mpi的基本实现

四.使用mpi的cell实现

五.使用mpi的异构实现

六.使用mpi的异构适用性实现

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include "mpi.h"
#include "malloc_align.h"

void environment_init();
typedef double data_t;
int N,M;//N行，M列，M = N+1
data_t * matrix_A;//X的系数矩阵
data_t * B;//y向量
data_t * X;//x的解向量
int *shift_table;
MPI_Status status;

//通过A(i,j)得到，x系数矩阵中[i,j]元素的地址
data_t * A(int i,int j){
    return matrix_A+i*M+j;
}
void print_matrix();

int main(int argc,char **argv)
{
    FILE * f;
    int i = 0,j = 0,index = 0,row = 0;
    int group_size;
    int *row_rank_table;
    int process_rownum;
    int myrank;
    data_t * mainrow_space;
    data_t ** row_address_table;


    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&group_size);
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
    //printf("mpi:%d %d\n",group_size,myrank);

    if(myrank == 0){
       f= fopen("/home/duanple/CBLAS/data.txt","r");
       fscanf(f,"%d%d",&N,&M);
       matrix_A = (data_t *)_malloc_align(sizeof(data_t)*N*M,7);
       B = (data_t *)_malloc_align(sizeof(data_t)*N,7);
       X = (data_t *)_malloc_align(sizeof(data_t)*N,7);
       //get matrix from file
       for(i = 0 ; i < N ; i++){
             for(j = 0 ; j < N ; j++){
               fscanf(f,"%lf",A(i,j));
        }
        fscanf(f,"%lf",A(i,j));
   *(B+i) = *A(i,j);
       }
   //print_matrix();
    }
    //broadcast the process row to be in,and rows num in every process
    MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
    M = N+1;
    //printf("rank:%d m:%d\n",myrank,M);
    //malloc the memory for the row of matrix
    shift_table = (int*)_malloc_align(sizeof(int)*N,7);
    row_rank_table = (int *)_malloc_align(sizeof(int)*N,7);
    row_address_table = (data_t**)_malloc_align(sizeof(data_t*)*N,7);
    //init shift table
    for(i = 0 ; i < N ; i++){
        shift_table[i] = i;
    }

    //divide matrix
    if(myrank == 0){
   for(i = 0,j = 0 ; i < N ; i++,j++){
        row_rank_table[i] = j%group_size;
   }
    }
    MPI_Bcast(row_rank_table,N,MPI_INT,0,MPI_COMM_WORLD);

    process_rownum = 0;
    for(i = 0 ; i < N ; i++){
   //printf("rank:%d row%d\n",myrank,row_rank_table[i]);
   if(row_rank_table[i] == myrank){
        row_address_table[i] = (data_t *)_malloc_align(sizeof(data_t)*M,7);
        process_rownum++;
   }
    }
    //
    //printf("rank:%d rows:%d\n",myrank,process_rownum);
    //root send the data to the processes.
    if(myrank == 0){//send
        for(i = 0 ; i < N ; i++){
       if(row_rank_table[i] == 0){
       row_address_table[i] = A(i,0);
       }else{
            MPI_Send(A(i,0),N,MPI_DOUBLE,row_rank_table[i],i,MPI_COMM_WORLD);
            MPI_Send(B+i,1,MPI_DOUBLE,row_rank_table[i],i,MPI_COMM_WORLD);
       }
        }
    }else{//receive
        for(i = 0 ; i < N ; i++){
       if(myrank !=0 && myrank == row_rank_table[i]){
           MPI_Recv(row_address_table[i],N,MPI_DOUBLE,0,i,MPI_COMM_WORLD,&status);
           MPI_Recv(row_address_table[i]+N,1,MPI_DOUBLE,0,i,MPI_COMM_WORLD,&status);
       }
        }
    }

   //printf the revecive result.
   /*
   for(i = 0 ; i < N ; i++){
   if(myrank != row_rank_table[i]) continue;
    printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
   for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
   }
   printf("\n");
    }   */
    //computer
    //消元
    mainrow_space = (data_t *)_malloc_align(sizeof(data_t)*M,7);//存放主行的临时空间
    for(i = 0 ; i < N-1 ; i++){
   if(myrank == row_rank_table[i]){

//上面的条件，可能造成同一个进程上的其他行不被处理。

           //选主元
       //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
       index = cblas_idamax(N-i,row_address_table[i]+i,1)+i;
           //列交换
                if(index != i){
                        int temp = shift_table[i];
              shift_table[i] = shift_table[index];
                  shift_table[index] = temp;

           data_t data_temp = *(row_address_table[i]+i);
           *(row_address_table[i]+i) = *(row_address_table[i]+index);
           *(row_address_table[i]+index) = data_temp;
           }
       //拷贝主行元素到临时空间
       //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
       memcpy(mainrow_space,row_address_table[i],sizeof(data_t)*M);
       //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
   }
   //其他行的列交换
       MPI_Bcast(&index,1,MPI_INT,row_rank_table[i],MPI_COMM_WORLD);
   if(myrank != row_rank_table[i]){

//上面的条件，可能造成同一个进程上的其他行不被处理。
                if(index != i){
                        int temp = shift_table[i];
              shift_table[i] = shift_table[index];
                  shift_table[index] = temp;
           //swap the other rows
           for(j = 0 ; j < N ; j++){
                if(myrank == row_rank_table[j]){
               data_t data_temp = *(row_address_table[j]+i);
               *(row_address_table[j]+i) = *(row_address_table[j]+index);
               *(row_address_table[j]+index) = data_temp;
                }
           }//end of for
           }

   }

   //发送主行元素
       MPI_Bcast(mainrow_space,M,MPI_DOUBLE,row_rank_table[i],MPI_COMM_WORLD);
        //主元归一化
   if(myrank == row_rank_table[i]){

//注意，上面这个条件是错的，因为当前rank上还有其他行，这样会导致其他行不被处理。
        //if(*(row_address_table[i]+i) == 0){continue;}
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",mainrow_space[j]);}printf("\n");

   }else{
        for(j = 0 ; j < N ; j++){
       //处理该进程所分配的行。
       if(j != i && myrank == row_rank_table[j]){
           if(*(row_address_table[j]+i) == 0){continue;}
              // 消元
               data_t factor = -1*(*(row_address_table[j]+i)/mainrow_space[i]);
           cblas_daxpy(N,factor,mainrow_space,1,row_address_table[j],1);
           *(row_address_table[j]+N) = *(row_address_table[j]+N) +mainrow_space[N]*factor;//elem of B
       }

        }

   }
    }

   //

   //printf the revecive result.
   /*
   for(i = 0 ; i < N ; i++){
   if(myrank != row_rank_table[i]) continue;
    printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
   for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
   }
   printf("\n");
    }   */

    //回代
    data_t sum = 0;

//注意一个进程不只一个行，因此只安排一个sum也是不够的。
    for(i = N-1 ; i>= 0 ; i--){
   if(myrank == row_rank_table[i]){
        mainrow_space[i] = (*(row_address_table[i]+N)-sum)/(*(row_address_table[i]+i));//(*(B+i) - sum) / *A(i,i);
   }
       MPI_Bcast(mainrow_space+i,1,MPI_DOUBLE,row_rank_table[i],MPI_COMM_WORLD);
   for(j = 0 ; j < i ; j++){
   //处理该进程所分配的行。
        if(myrank == row_rank_table[j]){
          sum += *(row_address_table[j]+i)*mainrow_space[i];
        }

   }


    }
    //printf res
    for(i = 0 ; myrank == 0 && i < N ; i++){
        printf("x%d:%lf ",shift_table[i],mainrow_space[i]);
    }

    /*
    //computer
    for(i = 0 ; i < N-1 ; i++){
        //选主元
   index = cblas_idamax(N-i,A(i,i),1)+i;//cblas_idamax 选择最大绝对值的blas库
   //printf("index:%d\n",index);
        //列交换
        if(index != i){
           temp = shift_table[i];
       shift_table[i] = shift_table[index];
           shift_table[index] = temp;
           for(row = 0 ; row < N ; row++){
               data_t temp = *A(row,i);
               *(A(row,i)) = *A(row,index);
               *A(row,index) = temp;
           }
        }
   for(row = 0 ; row < N ; row++)
       printf("%d ",shift_table[row]);
   printf("\n");
        //主元归一化
        // 消元
        for(row = i+1 ; row < N ; row++){
        if(*A(row,i) == 0) continue;
        data_t factor = -1*(*A(row,i)/ *A(i,i));
        cblas_daxpy(N-i,factor,A(i,i),1,A(row,i),1);//cblas_daxpy( N, alpha,X,incX, Y, incY); Y=Y+alpha*X
        *(B+row) = *(B+row) + *(B+i)*factor;
        }
        //print_matrix();

    }
    //回代，i代表了Xi,col代表了列号，用来计算前面已经计算的xi值的sum
    for(i = N-1 ; i >= 0 ; i--){
        data_t sum = 0;
        for( j = i+1 ; j < N ; j++){
            sum += *A(i,j)*X[j];
        }
        X[i] = (*(B+i) - sum) / *A(i,i);
    }
    //恢复,其中shift table保存了原来的列号
    for(i = 0 ; i < N ; i++){
        printf("x%d:%lf ",shift_table[i],X[i]);
    }
               */
    MPI_Finalize();

}
void print_matrix(){
    int i,j;
    printf("matrix:\n");
    for(i = 0 ; i < N ; i++){
        for(j = 0 ; j < N ; j++){
                printf("%lf ",*A(i,j));
        }
        printf("%lf\n",*(B+i));
    }

}

修改后的实现：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include "mpi.h"
#include "malloc_align.h"

void environment_init();
typedef double data_t;
int N,M;//N行，M列，M = N+1
data_t * matrix_A;//X的系数矩阵
data_t * B;//y向量
data_t * X;//x的解向量
int *shift_table;
MPI_Status status;

//通过A(i,j)得到，x系数矩阵中[i,j]元素的地址
data_t * A(int i,int j){
    return matrix_A+i*M+j;
}
void print_matrix();

int main(int argc,char **argv)
{
    FILE * f;
    int i = 0,j = 0,index = 0,row = 0;
    int group_size;
    int *row_rank_table;
    int process_rownum;
    int myrank;
    data_t * mainrow_space;
    data_t ** row_address_table;


    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&group_size);
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
    //printf("mpi:%d %d\n",group_size,myrank);

    if(myrank == 0){
        f= fopen("/home/duanple/CBLAS/data.txt","r");
        fscanf(f,"%d%d",&N,&M);
        matrix_A = (data_t *)_malloc_align(sizeof(data_t)*N*M,7);
        B = (data_t *)_malloc_align(sizeof(data_t)*N,7);
        X = (data_t *)_malloc_align(sizeof(data_t)*N,7);
        //get matrix from file
        for(i = 0 ; i < N ; i++){
             for(j = 0 ; j < N ; j++){
                fscanf(f,"%lf",A(i,j));
        }
        fscanf(f,"%lf",A(i,j));//A(i,N) also stores B[i]
    *(B+i) = *A(i,j);
        }
    //print_matrix();
    }
    //broadcast the process row to be in,and rows num in every process
    MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
    M = N+1;
    //printf("rank:%d m:%d\n",myrank,M);
    //malloc the memory for the row of matrix
    shift_table = (int*)_malloc_align(sizeof(int)*N,7);
    row_rank_table = (int *)_malloc_align(sizeof(int)*N,7);
    row_address_table = (data_t**)_malloc_align(sizeof(data_t*)*N,7);
    //init shift table
    for(i = 0 ; i < N ; i++){
        shift_table[i] = i;
    }

    //divide matrix
    if(myrank == 0){
    for(i = 0,j = 0 ; i < N ; i++,j++){
        row_rank_table[i] = j%group_size;
    }
    }
    MPI_Bcast(row_rank_table,N,MPI_INT,0,MPI_COMM_WORLD);

    process_rownum = 0;
    for(i = 0 ; i < N ; i++){
    //printf("rank:%d row%d\n",myrank,row_rank_table[i]);
    if(row_rank_table[i] == myrank){
        row_address_table[i] = (data_t *)_malloc_align(sizeof(data_t)*M,7);
        process_rownum++;
    }
    }
    //
    //printf("rank:%d rows:%d\n",myrank,process_rownum);
    //root send the data to the processes.
    if(myrank == 0){//send
        for(i = 0 ; i < N ; i++){
       if(row_rank_table[i] == 0){
        row_address_table[i] = A(i,0);
       }else{
            MPI_Send(A(i,0),N,MPI_DOUBLE,row_rank_table[i],i,MPI_COMM_WORLD);
            MPI_Send(B+i,1,MPI_DOUBLE,row_rank_table[i],i,MPI_COMM_WORLD);
       }
        }
    }else{//receive
        for(i = 0 ; i < N ; i++){
       if(myrank !=0 && myrank == row_rank_table[i]){
           MPI_Recv(row_address_table[i],N,MPI_DOUBLE,0,i,MPI_COMM_WORLD,&status);
           MPI_Recv(row_address_table[i]+N,1,MPI_DOUBLE,0,i,MPI_COMM_WORLD,&status);
       }
        }
    }

   //printf the revecive result.
   /*
   for(i = 0 ; i < N ; i++){
    if(myrank != row_rank_table[i]) continue;
       printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
    for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
    }
    printf("\n");
    }    */
    //computer
    //消元
    mainrow_space = (data_t *)_malloc_align(sizeof(data_t)*M,7);//存放主行的临时空间
    for(i = 0 ; i < N-1 ; i++){
    if(myrank == row_rank_table[i]){
            //选主元
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
        index = cblas_idamax(N-i,row_address_table[i]+i,1)+i;
            //列交换
                if(index != i){
                        int temp = shift_table[i];
               shift_table[i] = shift_table[index];
                   shift_table[index] = temp;

            data_t data_temp = *(row_address_table[i]+i);
            *(row_address_table[i]+i) = *(row_address_table[i]+index);
            *(row_address_table[i]+index) = data_temp;
            }
        //拷贝主行元素到临时空间
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
        memcpy(mainrow_space,row_address_table[i],sizeof(data_t)*M);
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
    }
    //其他行的列交换
        MPI_Bcast(&index,1,MPI_INT,row_rank_table[i],MPI_COMM_WORLD);
    //修改列交换标志，以进程为单位，不是以行
    if(myrank != row_rank_table[i]){
                  int temp = shift_table[i];
           shift_table[i] = shift_table[index];
               shift_table[index] = temp;
    }
    for(j = 0 ; j < N ; j++){
        if(j != i && myrank == row_rank_table[j]){
            data_t data_temp = *(row_address_table[j]+i);
            *(row_address_table[j]+i) = *(row_address_table[j]+index);
            *(row_address_table[j]+index) = data_temp;
        }
    }

    //发送主行元素
        MPI_Bcast(mainrow_space,M,MPI_DOUBLE,row_rank_table[i],MPI_COMM_WORLD);
    //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",mainrow_space[j]);}printf("\n");
        //主元归一化
    for(j = i+1 ; j < N ; j++){
    //处理该进程所分配的行。
        if(myrank == row_rank_table[j]){
            if(*(row_address_table[j]+i) == 0){continue;}
               // 消元
            data_t factor = -1*(*(row_address_table[j]+i)/mainrow_space[i]);
            cblas_daxpy(N,factor,mainrow_space,1,row_address_table[j],1);
            *(row_address_table[j]+N) = *(row_address_table[j]+N) +mainrow_space[N]*factor;//elem of B
        }

    }


    }

   //

   //printf the result.
   /*
   for(i = 0 ; i < N ; i++){
    if(myrank != row_rank_table[i]) continue;
       printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
    for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
    }
    printf("\n");
    }
    */
    //回代
    for(i = N-1 ; i>= 0 ; i--){
    if(myrank == row_rank_table[i]){
        mainrow_space[i] = *(row_address_table[i]+N)/(*(row_address_table[i]+i));//(*(B+i) - sum) / *A(i,i);
    }
        MPI_Bcast(mainrow_space+i,1,MPI_DOUBLE,row_rank_table[i],MPI_COMM_WORLD);
    for(j = 0 ; j < i ; j++){
    //处理该进程所分配的行。
        if(myrank == row_rank_table[j]){
          *(row_address_table[j]+N) -= *(row_address_table[j]+i)*mainrow_space[i];
        }

    }


    }
    //printf res

    for(i = 0 ; myrank == 0 && i < N ; i++){
        printf("x%d:%lf ",shift_table[i],mainrow_space[i]);
    }
    MPI_Finalize();

}
void print_matrix(){
    int i,j;
    printf("matrix:\n");
    for(i = 0 ; i < N ; i++){
        for(j = 0 ; j < N ; j++){
                printf("%lf ",*A(i,j));
        }
        printf("%lf\n",*(B+i));
    }

}

#############################################################################

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include "mpi.h"
#include "malloc_align.h"

void environment_init();
typedef double data_t;
int N,M;//N行，M列，M = N+1
data_t * matrix_A;//X的系数矩阵
data_t * B;//y向量
data_t * X;//x的解向量
int *shift_table;
MPI_Status status;

//通过A(i,j)得到，x系数矩阵中[i,j]元素的地址
data_t * A(int i,int j){
    return matrix_A+i*M+j;
}
void print_matrix();
double * processor_test(MPI_Comm comm,void test(),long loops,double comp_size,int master_node);
void data_malloc(double * processor_descrip,int * rank_table,int matrix_size);
void linear_solver(MPI_Comm comm);

void test(){
    double a = 1.0;
    double b = 2.44;
    double c = a*b;

}
int main(int argc,char **argv)
{


    MPI_Init(&argc,&argv);
    processor_test(MPI_COMM_WORLD,test,1000,1.0,0);
    //linear_solver(MPI_COMM_WORLD);
    MPI_Finalize();

}

//test the computer ability of processors
double * processor_test(MPI_Comm comm,void test(),long loops,double comp_size,int master_node){
    int myrank;
    int group_size,i;
    double start_time,end_time,time;
    double * processor_able;
    MPI_Comm_size(comm,&group_size);
    MPI_Comm_rank(comm,&myrank);
    if(myrank == master_node){
        processor_able = (double *)_malloc_align(sizeof(double)*group_size+2,7);
        processor_able[0] = group_size;
        processor_able[1] = comp_size*loops;
    } else
        processor_able = NULL;

    start_time = MPI_Wtime();
    for(i = 0; i < loops ; i++)
        test();
    end_time = MPI_Wtime();
    time = end_time-start_time;
    printf("rank%d %lf\n",myrank,time);

    MPI_Gather(&time,1,MPI_DOUBLE,processor_able+2,1,MPI_DOUBLE,master_node,comm);

    //print the test result.
    if(myrank == master_node){
        for(i = 0 ; i < group_size ; i++){
            printf("%d %lf ",i,processor_able[i+2]);
        }

    }

    return processor_able;

}
void data_malloc(double * processor_descrip,int * rank_table,int matrix_size){
    int p_num =(int)processor_descrip[0];


}
void linear_solver(MPI_Comm comm){

    FILE * f;
    int i = 0,j = 0,index = 0,row = 0;
    int group_size;
    int *row_rank_table;
    int process_rownum;
    int myrank;
    data_t * mainrow_space;
    data_t ** row_address_table;

    MPI_Comm_size(comm,&group_size);
    MPI_Comm_rank(comm,&myrank);
    //printf("mpi:%d %d\n",group_size,myrank);

    if(myrank == 0){
        f= fopen("/home/duanple/CBLAS/data.txt","r");
        fscanf(f,"%d%d",&N,&M);
        matrix_A = (data_t *)_malloc_align(sizeof(data_t)*N*M,7);
        B = (data_t *)_malloc_align(sizeof(data_t)*N,7);
        X = (data_t *)_malloc_align(sizeof(data_t)*N,7);
        //get matrix from file
        for(i = 0 ; i < N ; i++){
             for(j = 0 ; j < N ; j++){
                fscanf(f,"%lf",A(i,j));
        }
        fscanf(f,"%lf",A(i,j));//A(i,N) also stores B[i]
    *(B+i) = *A(i,j);
        }
    //print_matrix();
    }
    //broadcast the process row to be in,and rows num in every process
    MPI_Bcast(&N,1,MPI_INT,0,comm);
    M = N+1;
    //printf("rank:%d m:%d\n",myrank,M);
    //malloc the memory for the row of matrix
    shift_table = (int*)_malloc_align(sizeof(int)*N,7);
    row_rank_table = (int *)_malloc_align(sizeof(int)*N,7);
    row_address_table = (data_t**)_malloc_align(sizeof(data_t*)*N,7);
    //init shift table
    for(i = 0 ; i < N ; i++){
        shift_table[i] = i;
    }

    //divide matrix
    if(myrank == 0){
    for(i = 0,j = 0 ; i < N ; i++,j++){
        row_rank_table[i] = j%group_size;
    }
    }
    MPI_Bcast(row_rank_table,N,MPI_INT,0,comm);

    process_rownum = 0;
    for(i = 0 ; i < N ; i++){
    //printf("rank:%d row%d\n",myrank,row_rank_table[i]);
    if(row_rank_table[i] == myrank){
        row_address_table[i] = (data_t *)_malloc_align(sizeof(data_t)*M,7);
        process_rownum++;
    }
    }
    //
    //printf("rank:%d rows:%d\n",myrank,process_rownum);
    //root send the data to the processes.
    if(myrank == 0){//send
        for(i = 0 ; i < N ; i++){
       if(row_rank_table[i] == 0){
        row_address_table[i] = A(i,0);
       }else{
            MPI_Send(A(i,0),N,MPI_DOUBLE,row_rank_table[i],i,comm);
            MPI_Send(B+i,1,MPI_DOUBLE,row_rank_table[i],i,comm);
       }
        }
    }else{//receive
        for(i = 0 ; i < N ; i++){
       if(myrank !=0 && myrank == row_rank_table[i]){
           MPI_Recv(row_address_table[i],N,MPI_DOUBLE,0,i,comm,&status);
           MPI_Recv(row_address_table[i]+N,1,MPI_DOUBLE,0,i,comm,&status);
       }
        }
    }

   //printf the revecive result.
   /*
   for(i = 0 ; i < N ; i++){
    if(myrank != row_rank_table[i]) continue;
       printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
    for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
    }
    printf("\n");
    }    */
    //computer
    //消元
    mainrow_space = (data_t *)_malloc_align(sizeof(data_t)*M,7);//存放主行的临时空间
    for(i = 0 ; i < N-1 ; i++){
    if(myrank == row_rank_table[i]){
            //选主元
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
        index = cblas_idamax(N-i,row_address_table[i]+i,1)+i;
            //列交换
                if(index != i){
                        int temp = shift_table[i];
               shift_table[i] = shift_table[index];
                   shift_table[index] = temp;

            data_t data_temp = *(row_address_table[i]+i);
            *(row_address_table[i]+i) = *(row_address_table[i]+index);
            *(row_address_table[i]+index) = data_temp;
            }
        //拷贝主行元素到临时空间
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
        memcpy(mainrow_space,row_address_table[i],sizeof(data_t)*M);
        //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",*(row_address_table[i]+j));}printf("\n");
    }
    //其他行的列交换
        MPI_Bcast(&index,1,MPI_INT,row_rank_table[i],comm);
    //修改列交换标志，以进程为单位，不是以行
    if(myrank != row_rank_table[i]){
                  int temp = shift_table[i];
           shift_table[i] = shift_table[index];
               shift_table[index] = temp;
    }
    for(j = 0 ; j < N ; j++){
        if(j != i && myrank == row_rank_table[j]){
            data_t data_temp = *(row_address_table[j]+i);
            *(row_address_table[j]+i) = *(row_address_table[j]+index);
            *(row_address_table[j]+index) = data_temp;
        }
    }

    //发送主行元素
        MPI_Bcast(mainrow_space,M,MPI_DOUBLE,row_rank_table[i],comm);
    //printf("myrank:%d\n",myrank);for(j = 0 ; j < M ; j++){printf("%lf ",mainrow_space[j]);}printf("\n");
        //主元归一化
    for(j = i+1 ; j < N ; j++){
    //处理该进程所分配的行。
        if(myrank == row_rank_table[j]){
            if(*(row_address_table[j]+i) == 0){continue;}
               // 消元
            data_t factor = -1*(*(row_address_table[j]+i)/mainrow_space[i]);
            cblas_daxpy(N,factor,mainrow_space,1,row_address_table[j],1);
            *(row_address_table[j]+N) = *(row_address_table[j]+N) +mainrow_space[N]*factor;//elem of B
        }

    }


    }

   //

   //printf the result.
   /*
   for(i = 0 ; i < N ; i++){
    if(myrank != row_rank_table[i]) continue;
       printf("rank%d-row%d:\n",myrank,row_rank_table[i]);
    for(j = 0 ; j < M ; j++){
       printf("%lf ",*(row_address_table[i]+j));
    }
    printf("\n");
    }
    */
    //回代
    for(i = N-1 ; i>= 0 ; i--){
    if(myrank == row_rank_table[i]){
        mainrow_space[i] = *(row_address_table[i]+N)/(*(row_address_table[i]+i));//(*(B+i) - sum) / *A(i,i);
    }
        MPI_Bcast(mainrow_space+i,1,MPI_DOUBLE,row_rank_table[i],comm);
    for(j = 0 ; j < i ; j++){
    //处理该进程所分配的行。
        if(myrank == row_rank_table[j]){
          *(row_address_table[j]+N) -= *(row_address_table[j]+i)*mainrow_space[i];
        }

    }


    }
    //printf res

    for(i = 0 ; myrank == 0 && i < N ; i++){
        printf("x%d:%lf ",shift_table[i],mainrow_space[i]);
    }

}
void print_matrix(){
    int i,j;
    printf("matrix:\n");
    for(i = 0 ; i < N ; i++){
        for(j = 0 ; j < N ; j++){
                printf("%lf ",*A(i,j));
        }
        printf("%lf\n",*(B+i));
    }

}

评论这张

转发至微博

阅读(1109)| 评论(0)

历史上的今天

this.p={  m:2,
              b:2,
              loftPermalink:'',
              id:'fks_082066080084082068081086080095083086088068085081080068',
              blogTitle:'线性方程组求解-高斯消元',
              blogAbstract:'<p style=\"font-weight: bold; color: rgb(0, 128, 0);\"\><br\>4    5<br\>1.000000    4.000000    -2.000000    3.000000    6.000000<br\>2.000000    2.000000    0.000000    4.000000    2.000000<br\>3.000000    0.000000    -1.000000    2.000000    1.000000<br\>1.000000    2.000000    2.000000    -3.000000    8.000000<br\></p\>',
              blogTag:'',
              blogUrl:'blog/static/70971767200932294543260',
              isPublished:1,
              istop:false,
              type:2,
              modifyTime:1259251171000,
              publishTime:1240407943260,
              permalink:'blog/static/70971767200932294543260',
              commentCount:0,
              mainCommentCount:0,
              recommendCount:0,
              bsrk:-100,
              publisherId:0,
              recomBlogHome:false,
              currentRecomBlog:false,
              attachmentsFileIds:[],
              vote:{},
              groupInfo:{},
              friendstatus:'none',
              followstatus:'unFollow',
              pubSucc:'',
              visitorProvince:'',
              visitorCity:'',
              visitorNewUser:false,
              postAddInfo:{},
              mset:'000',
              mcon:'',
              srk:-100,
              remindgoodnightblog:false,
              isBlackVisitor:false,
              isShowYodaoAd:false,
              hostIntro:'',
              hmcon:'1',
              selfRecomBlogCount:'0',
              lofter_single:'<iframe width="140" height="560" style="overflow:hidden;" src="http://www.lofter.com/mailEntry.do?blogad=1&blog" frameBorder="0"></iframe>'
            }

{list a as x}
    {if !!x}
    <div class="iblock nbw-fce nbw-f40">
      <a class="fc03 noul" target="_blank" hidefocus="true" href="http://blog.163.com/${x.visitorName}/">
      {if x.visitorName==visitor.userName}
      <img alt="${x.visitorNickname|escape}" onerror="this.src=location.f40" class="cwd bdwa bdc0" src="${fn1(x.visitorName)}&r=${visitor.imageUpdateTime}"/>
      {else}
      <img alt="${x.visitorNickname|escape}" onerror="this.src=location.f40" class="cwd bdwa bdc0" src="${fn1(x.visitorName)}"/>
      {/if}
      </a>
      <div class="cwd vname thide">
        {if x.moveFrom=='wap'}
          <a class="noul pnt" target="_blank" href="http://blog.163.com/services/wapblog.html?frompersonalbloghome"><span title="来自网易手机博客" class="iblock wapIcon"> </span></a>
        {elseif x.moveFrom=='iphone'}
          <a class="noul pnt" target="_blank"><span title="来自iPhone客户端" class="iblock iphoneIcon"> </span></a>
        {elseif x.moveFrom=='android'}
          <a class="noul pnt" target="_blank"><span title="来自Android客户端" class="iblock androidIcon"> </span></a>
        {elseif x.moveFrom=='mobile'}
          <a class="noul pnt" target="_blank" href="http://blog.163.com/services/emsblog.html?frompersonalbloghome"><span title="来自网易短信写博" class="iblock wapIcon"> </span></a>
        {/if}
        <a class="fc03 m2a"  target="_blank" hidefocus="true" href="http://blog.163.com/${x.visitorName}/">
          ${fn(x.visitorNickname,8)|escape}
        </a>
      </div>
    </div>
    {/if}
    {/list}

<#--最新日志，群博日志--> <#--推荐日志-->

<p class="fc06">推荐过这篇日志的人：</p>
    <div>
      {list a as x}
      {if !!x}
      <div class="iblock nbw-fce nbw-f40">
        <a class="fc03 noul" target="_blank" hidefocus="true" href="http://blog.163.com/${x.recommenderName}/">
        <img alt="${x.recommenderNickname|escape}" onerror="this.src=location.f40" class="cwd bdwa bdc0" src="${fn1(x.recommenderName)}"/>
        </a>
        <div class="cwd thide">
          <a class="fc03 m2a" target="_blank" hidefocus="true" href="http://blog.163.com/${x.recommenderName}/">
            ${fn(x.recommenderNickname,6)|escape}
          </a>
        </div>
      </div>
      {/if}
      {/list}
    </div>
    {if !!b&&b.length>0}
    <p  class="fc06">他们还推荐了：</p>
    <ul>
    {list b as y}
      {if !!y}
        <li class="rrb"><span class="iblock">·</span><a class="fc03 m2a" target="_blank" href="http://blog.163.com/${y.recommendBlogPermalink}/?from=blog/static/70971767200932294543260">${y.recommendBlogTitle|escape}</a></li>
      {/if}
    {/list}
    </ul>
    {/if}

<#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇，下一篇--> <#-- 热度 -->

{list a as x}
    {if !!x}
    <div class="hotItem iblock nbw-fce nbw-f40">
      <a class="fc03 noul" target="_blank" hidefocus="true" href="http://blog.163.com/${x.publisherUsername}/">
      {if x.publisherUsername==visitor.userName}
      <img alt="${x.publisherNickname|escape}" onerror="this.src=location.f40" class="cwd bdwa bdc0" src="${fn1(x.publisherUsername)}&r=${visitor.imageUpdateTime}"/>
      {else}
      <img alt="${x.publisherNickname|escape}" onerror="this.src=location.f40" class="cwd bdwa bdc0" src="${fn1(x.publisherUsername)}"/>
      {/if}
      </a>
      <div class="cwd vname thide">
        <a class="fc03 m2a"  target="_blank" hidefocus="true" href="http://blog.163.com/${x.publisherUsername}/">
          ${fn(x.publisherNickname,8)|escape}
        </a>
      </div>
      <a class="f-myLikeIcons hottype {if x.type==1} js-liketype{elseif x.type==2} js-reblogtype{elseif x.type==3} js-sharetype{else}{/if}" target="_blank" hidefocus="true" href="http://blog.163.com/${x.publisherUsername}/"> </a>
    </div>
    {/if}
    {/list}

<#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->

页脚

我的照片书 - 手机博客 - 下载LOFTER APP - 订阅此博客

银河里的星星

导航

日志

线性方程组求解-高斯消元

历史上的今天

最近读者

热度

评论

页脚