1. Matlab solution Ax=b
Only the simplest way is written, where b needs to be a column vector, and semicolons are used to separate elements;
Ax = b
where the matrix A(nxn) and the column vector b(n) are known; the solution of the system of equations is stored in x, written in 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 iterative algorithm example
2.1 Principle
The principle of the iterative method is to create a convergent sequence of n-dimensional infinite terms. The convergence limit point x satisfies the equation Ax = b;
There are some techniques for using matrix A to create a convergent sequence. In the example here, since A is strictly dominated by the pivot element, it is easier to create convergence and the diagonal elements are directly removed;
For detailed principles, please refer to the numerical calculation textbook. It uses the induced norm of the matrix, spectral radius and other related theorems, and it is not complicated;
In order to express the principle here, I used the algorithm implementation in matrix language. The general principle is as follows. The principle is not written in a standardized way, so you can just skip it:
/**********************************
*
* ******Ax = b**********
* x = Bx + f
* x = -D'(L + U)x + D’b
*
* Dx = b – (L + U)x
*************************The following is the iterative process
* X1=(L + U) * X0
* X1=b-X1
* Dx = X1 ***(D here is a diagonal matrix, so you can get x by directly doing scalar division)
* x = D^-1 * X1
*X0 = x
*
*******************************/
The implementation of matrix language has two advantages. One is to facilitate traditional theoretical analysis, such as convergence analysis, and the other is to facilitate GPU parallel implementation;
2.2 Source code implementation
Here is the simplest Jacobi iteration method to find the solution of Ax=b;
gcc can be directly compiled and run;
jacobi.cpp
Source code:
#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 =\\ "); print_matrix(LpU, N, N); set_D(D, A, N); printf("D =\\ "); print_vector(D, N); memset(X0, 0, N*sizeof(float)); printf("X0 =\\ "); 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("\\ "); for(int i=0; i<N; i + + ){ for(int j=0; j<N; j + + ){ printf("%8.3f ", A[i + j*lda]); } printf("\\ "); } printf("\\ "); } void print_vector(float *V, int N){ printf("\\ "); for(int i=0; i<N; i + + ){ printf("%f ", V[i]); } printf("\\ "); } //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("\\ A =\\ "); print_matrix(A, N, N); printf("\\ y =\\ "); print_vector(y, N); sgemv(1, N, a, A, lda, b, y); printf("\\ y=Ax + y=\\ "); 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 =\\ "); 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; }
Compile:
Makefile
Jacobi: hello_Jacobi.cpp g + + -g $< -o $@ .PHONY: clean clean: -rm -rf Jacobi
As a reminder, the g++ line is preceded by a tab space, which is the syntax requirement of the Makefile, and the -rm is also preceded;
run: