1.  matlab 解 Ax=b

只写了最简单的方式，其中b需要是列向量，用分号隔开元素；

Ax = b

其中矩阵A(nxn) 和列向量 b(n) 已知；方程组的解存储在 x中，在matlab中写作：

> x=A\b

octave:7> A=[1,2; 1.0001, 2;]
A =

   1.0000   2.0000
   1.0001   2.0000

octave:8> b=[3; 3.0001;]
b =

   3.0000
   3.0001

octave:9> x=A\b
x =

   1.0000
   1.0000

octave:10> b-A*x
ans =

   0
   0

octave:11>

2. Ax=b 迭代法算法示例

2.1 原理

迭代法的原理，是创建一个n维的无穷项的收敛数列，收敛的极限点 x 满足方程 Ax = b ；

使用矩阵A来制造收敛数列，有一些技巧，这里的示例由于 A 是严格主元占优的，所以比较容易制造收敛，直接挖掉了对角元；

详细原理可以参考数值计算教材，用到了矩阵的诱导范数、谱半径等相关的定理，并不复杂；

这里为了表述原理，我采用了矩阵语言的算法实现，大概原理如下，原理写的不太规范，凑活着跳过即可：

/********************************
 *
 *   ******Ax = b*********
 *   x = Bx + f
 *   x = -D'(L+U)x + D'b
 *
 *  Dx = b - (L+U)x
 *  ***********************下面是迭代过程
 *  X1=(L+U) * X0
 *  X1=b-X1
 *  Dx = X1    ***(这里D是对角阵，所以直接做标量除法即可得到x)

*   x = D^-1   *  X1
 *  X0 = x
 *
 * ******************************/

矩阵语言的实现有两个好处，一个是便于传统的理论分析，比如收敛性分析等，另一个好处是有利于 gpu并行实现；

2.2 源码实现

这里给出最简单的 Jacobi 迭代法求 Ax=b 的解；

gcc 直接编译即可运行；

jacobi.cpp

源码：

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#define NA 3

void test_sgemv(int op, int N, float a, float *AA, int lda, float b, float *yy);
void sgemv(int op, int N, float a, float *A, int lda, float b, float *y);       // y = a*op(A)*X + b*y
void print_matrix(float *A, int N, int lda);
void print_vector(float *V, int N);
void set_LpU(float *LpU, float *A, int N);
void set_D(float *D, float *A, int N);
void Saxpy(int n, float alpha, float *x, float beta, float *y); //y = a*x + b*y
void x_div_y(int N, float* x, float* y);//x[i] = x[i]/y[i];
void Scopy(int n, float* x, float* y);// x -> y



void Saxpy(int n, float alpha, float *x, float beta, float *y) //y = a*x + b*y
{
        for(int i=0; i<n; i++){
                y[i] = alpha*x[i] + beta*y[i];
        }
}

void x_div_y(int N, float* x, float* y){//x[i] = x[i]/y[i];
        for(int i=0; i<N; i++)
                x[i] = x[i]/y[i];
}

void Scopy(int n, float* x, float* y)// x -> y
{
        for(int i=0; i<n; i++)
                y[i] = x[i];
}


void Jacobi_Ax_b(float *A, float *b, float *X0, float *Xk, int N, float eps, int iter){

        float *LpU = nullptr;
        //float *X0  = nullptr;
        //float *X1  = nullptr;
        float *D   = nullptr;

        LpU = (float*)malloc( N*N*sizeof(float));
        //X0  = (float*)malloc( N*sizeof(float));
        //X1  = (float*)malloc( N*sizeof(float));
        D   = (float*)malloc( N*sizeof(float));

        //1.0 set LpU, set D, X0 = 0;
        set_LpU(LpU, A, N);             printf("LpU =\n");      print_matrix(LpU, N, N);
        set_D(D, A, N);                 printf("D =\n");        print_vector(D, N);
        memset(X0, 0, N*sizeof(float)); printf("X0 =\n");       print_vector(Xk, N);

        //2.0 loop
/* S1  X1=LU* x0
 * S2  X1=b-X1
 * S3  Dx = X1
 * S4  X1 = X1/D
 */
        for(int it=0; it<iter; it++){
                //void sgemv(int op, int N, float a, float *A, int lda, float b, float *y);
                // S1  X1=LU* X0
                sgemv(0, N, 1.0f, LpU, N, 0.0f, Xk);
                //void Saxpy(int n, float alpha, float *x, float beta, float *y); //y = a*x + b*y
                // S2  X1=b-X1
                Saxpy(N, 1.0f, b, -1.0f, Xk);
                //void x_div_y(int N, float* x, float* y){//x[i] = x[i]/y[i];
                // S3 DXk+1 = Xk
                x_div_y(N, Xk, D);
        }
        //Scopy(N, X1, xk);

}

void set_LpU(float *LpU, float *A, int N){
        for(int i=0; i<N; i++){
                for(int j=0; j<N; j++){
                        LpU[i + j*N] = A[i + j*N];
                        if(i == j)
                                LpU[i+i*N] = 0.0f;
                }
        }
}

void set_D(float *D, float *A, int N){
        for(int i=0; i<N; i++)
        D[i] = A[i + i*N];
}

void sgemv(int op, int N, float a, float *A, int lda, float b, float *y){
        // y = a*op(A)*X + b*y
        float* y0 = nullptr;

        y0 = (float*)malloc(N*sizeof(float));
        memcpy(y0, y, N*sizeof(float));

        for(int i=0; i<N; i++){
                float sigma = 0.0f;

                for(int j=0; j<N; j++){
                        sigma += A[i + j*lda]*y0[j];
                }
                y[i] = a*sigma + b*y[i];
        }

        free(y0);
}

void print_matrix(float *A, int N, int lda){
        printf("\n");
        for(int i=0; i<N; i++){
                for(int j=0; j<N; j++){
                        printf("%8.3f ", A[i + j*lda]);
                }
                printf("\n");
        }
        printf("\n");
}

void print_vector(float *V, int N){
        printf("\n");
        for(int i=0; i<N; i++){
                printf("%f ", V[i]);
        }
        printf("\n");
}


//void sgemv(int op, int N, float a, float *A, int lda, float b, float *y){
void test_sgemv(int op, int N, float a, float *AA, int lda, float b, float *yy){
        float *y = nullptr;
        float *A = nullptr;

        y = (float*)malloc(N*sizeof(float));
        A = (float*)malloc(N*N*sizeof(float));

        memcpy(y, yy, N*sizeof(float));
        memcpy(A, AA, N*N*sizeof(float));

        printf("\nA =\n");
        print_matrix(A, N, N);
        printf("\ny =\n");
        print_vector(y, N);

        sgemv(1, N, a, A, lda, b, y);

        printf("\n y=Ax+y=\n");
        print_vector(y, N);

        free(y);
        free(A);
}



int main(){

        float A[NA*NA] =
        {
                10, -1, -1, -1, 10, -1, -2, -2, 5
        //      10, 2, 1, 3, -10, 3, 1, 3, 10// column major
        /*
                10,   3,  1,
                 2, -10,  3,
                 1,   3, 10*/
        };

        float b[NA]={
                7.2, 8.3, 4.2
                //14, -5, 14
                };

// void test_sgemv(int op, int N, float a, float *AA, int lda, float b, float *yy){
        //test_sgemv(1, NA, 1.0f, A, NA, 1.0f, b);

        float* Ah = nullptr;
        float* bh = nullptr;
        int N = NA;
        Ah = (float*)malloc(N*N*sizeof(float));
        bh = (float*)malloc(N*sizeof(float));

        memcpy(Ah, A, N*N*sizeof(float));
        memcpy(bh, b, N*sizeof(float));
        printf("Ah =\n");
        print_matrix(Ah, N, N);

        float *x0 = nullptr;
        float *x1 = nullptr;

        x0 = (float*)malloc(N*sizeof(float));
        x1 = (float*)malloc(N*sizeof(float));

        float eps = 1.0e-7;
        int iter  = 15;

        Jacobi_Ax_b(Ah, bh, x0, x1, N, eps, iter);

        print_vector(x1, N);

        return 0;
}

编译：

Makefile

Jacobi: hello_Jacobi.cpp
        g++ -g $< -o $@


.PHONY: clean
clean:
        -rm -rf Jacobi

絮叨地提示一下，g++ 行的前面是一个 tab空格，Makefile的语法要求，-rm 前面也是；

运行：

3.3 matlab 对结果进行验证：