CN106681960A - Acceleration method for solution of linear equation set with Monte Carlo method based on shared memory - Google Patents

Abstract

The invention provides an acceleration method for solution of a linear equation set with a Monte Carlo method based on shared memory. The method is characterized in that a sequence to be randomized is stored in the shared memory and the storage is conducted through optimization. Hence, Bank Conflict will not be generated from access of different threads to the shared memory; the reading speed is increased; and the whole algorithm is accelerated.

Description

An Accelerated Method for Solving Linear Equations Based on Shared Memory Monte Carlo Method

technical field

The invention relates to the field of numerical calculation, and more specifically, relates to an acceleration method for solving linear equations based on a shared memory Monte Carlo method.

Background technique

For large-scale scientific and engineering computing, because many problems can be attributed to linear equations, the solution method of large-scale sparse linear equations is one of the core issues, which has important theoretical significance and practical application value. The Monte Carlo method is a numerical calculation method based on the theory of probability and statistics, which can solve many practical problems by using random numbers. Since the Monte Carlo method was proposed to solve linear equations in the 1950s, many improvements have been made and it has been widely used in financial engineering, computational physics, and other fields.

At present, there are mainly the following methods for linear equations: (1) direct method, which is a method to obtain the real solution of linear equations directly through certain arithmetic operations, such as Gaussian elimination method and Clem method; ( 2) The iterative method is a method that starts from an initial estimate and gradually approaches the real solution through an iterative method during the solution process, such as Jacobi iterative method and Gauss-Seidel iterative method; (3) Monte Carlo method , is a probability method developed based on the Jacobi method, which simulates a process by using random numbers, and estimates the solution of the problem by finding the statistics of the mathematical expectation of its specific random variable.

At present, the calculation for solving large-scale sparse linear equations is huge, and it is not actually necessary to solve all the solutions. The Monte Carlo method can be used to solve a single approximate solution of large-scale sparse linear equations, and the convergence speed of the Monte Carlo method has nothing to do with the size of the linear equations, which is more efficient. At the same time, because the problem itself has natural parallelism, GPU-based parallel algorithms can be used to accelerate the solution of linear equations.

Contents of the invention

The invention provides an acceleration method for solving linear equations based on a shared memory Monte Carlo method, which can accelerate the speed of parallel algorithms and achieve higher efficiency.

In order to achieve the above-mentioned technical effect, the technical scheme of the present invention is as follows:

An acceleration method for solving a linear equation system based on a shared memory Monte Carlo method, comprising the following steps:

S1: Preprocess the input data to obtain the data needed to solve the linear equation system;

S2: According to the data obtained in S1, use the quasi-random number generation algorithm that comes with the cuRAND library in the CUDA architecture to generate simulation results that require quasi-random sequences and target solution components;

S3: Store the simulation results of the quasi-random sequence and the target solution components in the shared memory allocated by the GPU, use the GPU to solve them in parallel, and finally use the GPU to perform reduction and summation to obtain the required components.

Further, the process of step S1 is as follows:

S11: Input coefficient matrix A, vector b, truncation error δ, probability error ε, and the number of quasi-random sequences N is the number of Markov chains,

Among them, a linear equation system with n unknowns that needs to be solved is expressed as: Ax=b, A is a non-singular sparse matrix of n×n, b is a given known vector of n×1 dimension, and x is An n×1-dimensional unknown vector, which is the vector to be solved,

The following Jacobian iterative form is obtained by transformation: x=Gx+f, G is a matrix of n×n, and f is a vector of n×1;

S12: Calculate moment G:

Among them, a _ij is the element in the coefficient matrix A;

S13: Calculate the vector f:

Among them, a _ij is the element in the coefficient matrix A, and b _i is the element in the vector b;

S14: Calculate the row norm:

Among them, g _ij is an element in the matrix G;

S15: Calculate transition probability matrix P:

Further, the process of step S2 is as follows:

S21: Declare a random number pointer, and use the cuRAND library of the CUDA architecture in the GPU to specify a quasi-random number generation algorithm;

S22: Set a random number seed for the quasi-random sequence generation algorithm;

S23: Specify the step size of the quasi-random sequence, generate a corresponding number of quasi-random sequences, and save them in the global memory.

Further, the process of step S3 is as follows:

S31: first apply for the GPU to allocate shared memory, and store the quasi-random sequence generated in step S2 and a simulation result of the target solution component;

S32: use GPU to simulate each Markov chain in parallel;

S33: sum the Markov chain simulation results in all blocks and divide by the Markov chain number N, and finally obtain the value of the target component that needs to be solved, namely

Further, the process of step S32 is specifically as follows:

S321: For each thread warp, make each thread read the quasi-random sequence stored in the shared memory in parallel, and store a part of the quasi-random sequence in the shared memory corresponding to each block;

S322: Make the initial state the point corresponding to the component to be solved;

S323: Compare the corresponding quasi-random numbers in the quasi-random sequence with the state transition probability corresponding to the probability transition matrix, select the probability interval that the quasi-random number falls into, and select the next state corresponding to the transition probability as the Markov chain next state;

S324: simulate a Markov chain according to the C22 method, obtain the simulation result of a chain, and store it in the shared memory;

S325: sum the results of all Markov chains in each block and store in the shared memory.

Compared with the prior art, the beneficial effects of the technical solution of the present invention are:

The Monte Carlo method based on the shared memory proposed by the present invention is an acceleration method for solving linear equations. The quasi-random sequence is stored in the shared memory, and is stored through an optimization method, so that different threads do not have memory conflicts when accessing the shared memory. (Bank Conflict), which improves the reading speed and speeds up the entire algorithm.

Description of drawings

Fig. 1 is the schematic flow chart of the accelerated method for solving linear equations based on the Monte Carlo method of shared memory in the present invention;

Fig. 2 is the concrete storage mode of shared memory of the present invention;

Figure 3 is an illustration of shared memory bank conflicts.

detailed description

The accompanying drawings are for illustrative purposes only and cannot be construed as limiting the patent;

In order to better illustrate this embodiment, some parts in the drawings will be omitted, enlarged or reduced, and do not represent the size of the actual product;

For those skilled in the art, it is understandable that some well-known structures and descriptions thereof may be omitted in the drawings.

The technical solutions of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

Example 1

As shown in Figure 1, the present invention first wakes up the preprocessing of the input data, then uses the Kernel to start the GPU and allocates the corresponding required resources, then generates a quasi-random sequence with a fixed step size and number, and optimizes the quasi-random sequence into In the shared memory, use the thread warp to read the quasi-random sequence in the shared memory in parallel, and simulate the Markov chain to obtain its simulation results, and divide the sum of the simulation results of all Markov chains by the Markov chain The number of husband chains N to obtain the value of the target component that needs to be obtained

As shown in Figure 2, 32 threads form a thread warp, and the following is a shared storage, and the shared storage is composed of alternately arranged storage bodies (banks). Generally speaking, the number of threads in a warp is the same as the number of banks in shared storage, so threads in a warp can access data stored in shared storage in parallel. We store the quasi-random sequences with a step size of k into the same storage bank of the shared memory in sequence. Each thread in the same thread warp needs a complete quasi-random sequence to simulate a Markov chain, and a Threads only need to access different addresses of the same memory bank to complete a quasi-random sequence of reading. Different threads only need to access and read different memory banks. In GPU, different memory banks can simultaneously access and read Therefore, each thread in the same thread warp can optimize the stored quasi-random sequence from the present invention in parallel, that is, prepare for a parallel simulated Markov chain without bank conflict.

As shown in Figure 3, a shared memory consists of many alternately arranged memory banks. Assume that bank 0 (bank 0) in the figure, that is, the gray column, stores the first quasi-random number of a partial quasi-random sequence, so that different threads need to access different addresses of bank 0 (bank 0) at the same time , which will cause a bank conflict. However, in order to complete the reading of quasi-random numbers, multiple threads will be generated to request access to the same bank, but only one thread can access bank 0 (bank0) at a time, and multiple threads can only access serially to obtain results , resulting in a performance loss.

The specific steps of the acceleration method based on the Monte Carlo solution of the linear equation system based on the shared memory of the present invention are as follows:

1. First, input the data of the linear equation system to be calculated. Including the coefficient matrix A, the vector b, the truncation error δ of the Markov chain, the probability error ε, and the number N of the Markov chain, which is the number of quasi-random number sequences. At the same time, it is necessary to know the target component x _m to be solved Subscript Idx.

2. Data preprocessing part: In order to make the sequence of the original input form Ax=b change into the Jacobian form x=Gx+f through a point transformation, use the input coefficient matrix A and vector b to calculate the Jacobian form The matrix G and vector f of , the formulas are: Among them, a _ij is the element in the coefficient matrix A, and b _i is the element in the vector b, that is, the data in the Jacobian form is obtained through the data in the original form. Then the row norm of the matrix G can be calculated: Where g _ij is an element in the matrix G. After obtaining the row norm, you can calculate the probability transition matrix P that needs to be used later: Where g _is the element in the matrix G.

3. Generate a quasi-random number sequence. First declare the random number pointer, use the cuRAND library of the CUDA (Compute Unified Device Architecture, Unified Computing Architecture) architecture in the GPU (an algorithm library for generating your random number sequence), and specify the quasi-random number generation algorithm as needed. The random sequence generation algorithm sets the random number seed, and at the same time specifies the quasi-random sequence step size k, uses the quasi-random number generation algorithm to generate a quasi-random sequence with a number of N, and saves it in the global memory.

4. In order to achieve parallelization, it is necessary to use the Kernel to start the GPU and allocate resources to the corresponding blocks, such as the shared memory owned by each block. There are many threads in each block block, and every 32 threads form a thread warp, and work in units of thread warp.

5. Optimize the storage of quasi-random sequences into shared memory. First, the quasi-random sequence is read from the global memory, and many quasi-random sequences form a 2D data. As shown in Figure 2, the shared memory is composed of many memory banks (banks) arranged alternately. The quasi-random sequence and the same sequence are stored in the memory bank in columns, so that for a group of threads in a warp warp In other words, for example, thread 0 can read the quasi-random sequence with bank index 0 without restriction, without generating bank conflicts, so that a thread warp can simultaneously access multiple different memory banks to simultaneously read out the required quasi-random sequence without bank conflict, which can improve efficiency and obtain better parallelism.

6. After optimizing the shared memory in step 5, the threads in the warp can access and read quasi-random sequences in parallel. For each warp, each thread in it can read and store in parallel in the shared memory. The quasi-random sequence in the memory stores part of the quasi-random sequence in the shared memory corresponding to each block. After obtaining the quasi-random sequence, the Markov chain can be simulated. By obtaining the subscript Idx of the target component that needs to be solved, the initial state of the Markov chain is shared with the corresponding state (state←Idx) for the target to be solved, and the initial expectation E is 0, and then according to the quasi-random sequence In order, compare with the probability transition matrix P[state][j], j=1,...,n in order to obtain the probability interval that the quasi-random number falls into, and the corresponding state s is the Markov The next state of the chain (new_state←s), keeps looping until the simulation of a Markov chain is completed, and the simulation result is obtained. At the same time, the results of all Markov chains in each block are summed and stored in in shared memory.

7. Finally, after all (N) Markov chains are simulated, all the results are summed to get the average value, and finally the target component x _m to be solved can be obtained.

The same or similar reference numerals correspond to the same or similar components;

The positional relationship described in the drawings is only for illustrative purposes and cannot be construed as a limitation to this patent;

Apparently, the above-mentioned embodiments of the present invention are only examples for clearly illustrating the present invention, rather than limiting the implementation of the present invention. For those of ordinary skill in the art, other changes or changes in different forms can be made on the basis of the above description. It is not necessary and impossible to exhaustively list all the implementation manners here. All modifications, equivalent replacements and improvements made within the spirit and principles of the present invention shall be included within the protection scope of the claims of the present invention.

Claims (5)

1. an accelerated method for solving linear equations based on the Monte Carlo method of shared memory, it is characterized in that, comprising the following steps: S1: Preprocess the input data to obtain the data needed to solve the linear equation system; S2: According to the data obtained in S1, use the quasi-random number generation algorithm that comes with the cuRAND library in the CUDA architecture to generate simulation results that require quasi-random sequences and target solution components; S3: Store the simulation results of the quasi-random sequence and the target solution components in the shared memory allocated by the GPU, use the GPU to solve them in parallel, and finally use the GPU to perform reduction and summation to obtain the required components. 2. the method for accelerating the solving of linear equations based on the Monte Carlo method of shared memory according to claim 1, is characterized in that, the process of described step S1 is as follows: S11: Input coefficient matrix A, vector b, truncation error δ, probability error ε, and the number of quasi-random sequences N is the number of Markov chains, Among them, a linear equation system with n unknowns that needs to be solved is expressed as: Ax=b, A is a non-singular sparse matrix of n×n, b is a given known vector of n×1 dimension, and x is An n×1-dimensional unknown vector, which is the vector to be solved, The following Jacobian iterative form is obtained by transformation: x=Gx+f, G is a matrix of n×n, and f is a vector of n×1; S12: Calculate moment G: ${<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; g</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; \}</span>}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; =</span>\left\{\begin{array}{c}\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; w</span>}\mathrm{<span\; class="notranslate">\; h</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; j</span>}\\ \frac{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}}{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; w</span>}\mathrm{<span\; class="notranslate">\; h</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; \&NotEqual;</span>\mathrm{<span\; class="notranslate">\; j</span>}\end{array}\right.$ Among them, a _ij is the element in the coefficient matrix A; S13: Calculate the vector f: ${\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; no</span>}$ Among them, a _ij is the element in the coefficient matrix A, and b _i is the element in the vector b; S14: Calculate the row norm: $\mathrm{<span\; class="notranslate">\; R</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; w</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; g</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; no</span>}$ Among them, g _ij is an element in the matrix G; S15: Calculate transition probability matrix P: ${\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; g</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; |</span>}{\mathrm{<span\; class="notranslate">\; R</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; w</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; .</span>$ 3. the method for accelerating the solving of linear equations based on the Monte Carlo method of shared memory according to claim 2, is characterized in that, the process of described step S2 is as follows: S21: Declare a random number pointer, and use the cuRAND library of the CUDA architecture in the GPU to specify a quasi-random number generation algorithm; S22: Set a random number seed for the quasi-random sequence generation algorithm; S23: Specify the step size of the quasi-random sequence, generate a corresponding number of quasi-random sequences, and save them in the global memory. 4. the method for accelerating the solving of linear equations based on the Monte Carlo method of shared memory according to claim 3, is characterized in that, the process of described step S3 is as follows: S31: first apply for the GPU to allocate shared memory, and store the quasi-random sequence generated in step S2 and a simulation result of the target solution component; S32: use GPU to simulate each Markov chain in parallel; S33: sum the Markov chain simulation results in all blocks and divide by the Markov chain number N, and finally obtain the value of the target component that needs to be solved, namely 5. the method for accelerating the solution of linear equations based on the Monte Carlo method of shared memory according to claim 4, is characterized in that, the process of described step S32 is specifically as follows: S321: For each thread warp, make each thread read the quasi-random sequence stored in the shared memory in parallel, and store a part of the quasi-random sequence in the shared memory corresponding to each block; S322: Make the initial state the point corresponding to the component to be solved; S323: Compare the corresponding quasi-random numbers in the quasi-random sequence with the state transition probability corresponding to the probability transition matrix, select the probability interval that the quasi-random number falls into, and select the next state corresponding to the transition probability as the Markov chain next state; S324: simulate a Markov chain according to the C22 method, obtain the simulation result of a chain, and store it in the shared memory; S325: sum the results of all Markov chains in each block and store in the shared memory.