A heterogeneous scientific computing acceleration method based on a video memory zero-copy mechanism

Abstract

The application discloses a heterogeneous scientific calculation acceleration method based on a display memory zero-copy mechanism, and is applied to a heterogeneous calculation system equipped with a central processing unit and a graphic processing unit. In an initialization stage, a constraint matrix of a linear programming problem is transmitted to high-bandwidth display memory of the graphic processing unit in a compressed sparse row format at one time and a lock page residence operation is performed. In an iteration solving stage, the central processing unit only receives scalar-level summary information returned by the graphic processing unit to perform selection decisions of on-base variables and off-base variables, and sends byte-level control instructions to the graphic processing unit. The graphic processing unit autonomously completes all intensive operations of simple cost parallel calculation, minimum ratio parallel inspection and base matrix implicit updating operation in the display memory, so that the total bus transmission data amount is reduced from dozens of megabytes per iteration to dozens of bytes constant level, and the bandwidth bottleneck problem caused by frequent matrix data transfer in the existing heterogeneous calculation architecture is solved.

Description

Technical Field

[0001] This application relates to the fields of heterogeneous computing and numerical optimization technology, and in particular to a method for accelerating heterogeneous scientific computing based on a zero-copy mechanism of video memory. Background Technology

[0002] In heterogeneous computing architectures that utilize the CPU and GPU to collaboratively solve large-scale linear programming problems, each iteration of the modified simplex method requires transferring the column vector data of the basis matrix between the CPU's main memory and the GPU's video memory via an external component interconnect bus. Since the effective bandwidth of this bus is far lower than the access bandwidth of the GPU's high-bandwidth on-chip video memory, data transfer latency accounts for a major portion of the total time consumed in a single iteration, forming a bandwidth bottleneck that limits the solution throughput. Existing technologies alleviate this bottleneck by reducing the amount of data transmitted through sparse compression formats or by overlapping computation and communication pipelines to mask transmission latency, but both rely on the premise that matrix data must be transferred across processors in each iteration, failing to fundamentally eliminate the data transfer overhead. Summary of the Invention

[0003] Firstly, this application provides a heterogeneous scientific computing acceleration method based on a zero-copy mechanism for video memory, comprising: transferring the constraint matrix, right-hand constraint vector, and objective function cost vector of a linear programming problem from the CPU's main memory to the high-bandwidth video memory of the graphics processor in one go; performing a page-locking and resident operation on the video memory region storing the constraint matrix to make the constraint matrix resident during the iterative solution process; initializing the current basic feasible solution based on the right-hand constraint vector in the high-bandwidth video memory; performing a reduction cost operation on all non-basic variables based on the resident constraint matrix on the graphics processor to obtain the global minimum reduction cost value and the corresponding entering-basic column index, and transmitting it back to the CPU in the form of a scalar signal for judgment. Determine if the current solution satisfies the optimality condition. If it does not, extract the column vector corresponding to the in-basic column index at the graphics processor and perform a basis inverse product operation to obtain the updated column vector. Based on the updated column vector and the current basic feasible solution, perform a minimum ratio test to obtain the off-basic row index. The off-basic row index is then sent back to the central processing unit in the form of a scalar signal. Based on the in-basic column index and the off-basic row index, the basis matrix is ​​implicitly updated in-place in high-bandwidth video memory using a product form of basis inverse representation, and the basic feasible solution is updated. The central processing unit updates the basic variable index set and the non-basic variable index set and sends an instruction to the graphics processor to synchronize the copy of the basic variable index set. Then, it returns to perform a simplified cost operation until the optimality condition is satisfied.

[0004] Optionally, the one-time transfer uses a compressed sparse row format, including: performing a format conversion on the constraint matrix to generate a row pointer array, column index array, and value array in compressed sparse row format; calling the graphics processing unit device detection interface to obtain the total capacity of the high-bandwidth video memory, and estimating the required video memory based on the total number of non-zero elements, the number of constraint rows, and the number of decision variables in the constraint matrix; terminating execution when the required video memory exceeds a preset percentage threshold of the total capacity; allocating a page-locked memory region on the central processing unit and filling the row pointer array, the column index array, the value array, the right-hand constraint vector, and the objective function cost vector into the page-locked memory region; and writing all the data in the page-locked memory region to the pre-planned offset address of the high-bandwidth video memory through a one-time direct memory access transfer.

[0005] Optionally, the implicit update operation of the basis matrix includes: obtaining the pivot value of the updated column vector at the index position away from the basis row, and obtaining the reciprocal of the pivot value; performing a linear transformation on each row component of the updated column vector based on the reciprocal to generate an elementary transformation vector, and appending the elementary transformation vector to the transformation chain of the basis inverse representation.

[0006] Optionally, the basis inverse product operation uses the product form basis inverse representation, which expresses the inverse of the basis matrix as a chain product of multiple elementary transformation vectors arranged in order. The basis inverse product operation selects the corresponding solution process according to the type of the input vector: when the input is a column vector, a forward solution process is used, and a scattering rank-one update operation is performed on the input column vector in the order of generation of the elementary transformation vectors to obtain the product of the inverse of the basis matrix and the input column vector. Each step of the scattering rank-one update operation includes extracting the component value of the input column vector at the pivot position as... After scaling, the pivot component is cleared to zero, and the product of the scaling factor and the corresponding elementary transformation vector is added to all components of the input column vector. When the input is a row vector, a reverse solution process is adopted, and the inner product assignment operation is performed on the input row vector in reverse order of the generation order of the elementary transformation vectors to obtain the product of the input row vector and the inverse of the basis matrix. Each step of the inner product assignment operation includes calculating the inner product value of the input row vector and the current elementary transformation vector and assigning the inner product value to the component of the input row vector at the corresponding pivot position, while the other components of the input row vector remain unchanged.

[0007] Optionally, performing simplified cost calculation on all non-basic variables includes: extracting corresponding components from the objective function cost vector based on the basic variable index set through the reverse solution process to form the cost coefficient vector of the current basic variable; using the cost coefficient vector as the input row vector to obtain the dual variable vector; allocating a parallel execution thread for each non-basic variable; each parallel execution thread reading the non-zero elements of the corresponding constraint matrix column from the high-bandwidth video memory and performing a sparse inner product operation with the dual variable vector; and subtracting the result of the sparse inner product operation from the corresponding original cost coefficient to obtain the simplified cost value of the non-basic variable.

[0008] Optionally, the step of obtaining the global minimum parsimony cost value and the corresponding entering base column index adopts a two-level parallel reduction structure, including: in the first-level reduction, the parsimony cost values ​​of all non-basic variables obtained by the parsimony cost operation are evenly distributed to multiple thread blocks, and multiple threads in each thread block perform tree-structured comparison reduction in shared memory to obtain the local minimum value and its local index; in the second-level reduction, global reduction is performed on the local minimum values ​​output by all the thread blocks to obtain the global minimum parsimony cost value and the entering base column index.

[0009] Optionally, the extraction of the column vector corresponding to the infeasor column index includes: when the infeasor column index points to the original decision variable column, allocating a parallel execution thread to each row of the constraint matrix, and each parallel execution thread performing a binary search within the column index range of the compressed sparse row format of the corresponding row to locate the non-zero element value corresponding to the infeasor column index, and assembling the search results of all rows into a dense column vector; when the infeasor column index exceeds the range of the original decision variable column of the constraint matrix, the column vector corresponding to the infeasor column index is a standard basis vector, and the standard basis vector is directly generated as the dense column vector.

[0010] Optionally, the minimum ratio test adopts the Brand inverse loop rule. When there are multiple basic variable rows whose ratios are all equal to the minimum ratio, the row with the smallest index value in the basic variable index set is selected as the off-base row index.

[0011] Optionally, it also includes: monitoring the cumulative number of elementary transformation vectors in the product form basis inverse representation; when the cumulative number exceeds a preset reconstruction threshold, performing partial principal triangular decomposition on the current basis matrix to regenerate the product form basis inverse representation and resetting the cumulative number.

[0012] Optionally, it further includes: after each iteration of a preset verification interval, selecting a basic variable column index from the basic variable index set, performing the basis inverse product operation on the constraint matrix column vector corresponding to the basic variable column index through the forward solution process to obtain the verification vector, comparing the infinite norm residual between the verification vector and the standard basis vector corresponding to the position of the basic variable column index in the basis matrix with a preset precision threshold, and triggering the partial principal component triangular decomposition when the infinite norm residual exceeds the preset precision threshold.

[0013] Optionally, determining whether the current solution satisfies the optimality condition includes: comparing the global minimum simplified cost value with the negative numerical precision tolerance; when the global minimum simplified cost value is not less than the negative numerical precision tolerance, determining that the current solution satisfies the optimality condition.

[0014] Optionally, when the optimality condition is satisfied, the method further includes: traversing the basic variable index set to detect whether there are basic variables whose index values ​​exceed the range of the original decision variable columns of the constraint matrix and are still in the basic variable position; when there are basic variables that exceed the range and their corresponding basic feasible solution component values ​​are greater than the preset feasibility tolerance, it is determined that the linear programming problem has no feasible solution; when there are no basic variables that exceed the range or the basic feasible solution component values ​​corresponding to the basic variables that exceed the range are not greater than the preset feasibility tolerance, the current basic feasible solution is mapped to the complete decision variable solution vector according to the basic variable index set at the graphics processor end, and the complete decision variable solution vector is transferred back from the high-bandwidth video memory to the central processing unit main memory through a one-time direct memory access transfer.

[0015] Optionally, the scalar signal return is achieved by establishing an asynchronous direct memory access channel between the device memory on the graphics processor side and the page-locked memory on the central processing unit side, with each data return not exceeding 32 bytes.

[0016] Optionally, the page-locking residency operation includes calling the video memory suggestion interface to mark the video memory region as a read-only priority access mode, and calling the video memory prefetch interface to prefetch data pages in the video memory region into the graphics processor's L2 cache.

[0017] The beneficial effects of this application are as follows: This application loads the constraint matrix into the high-bandwidth video memory of the graphics processor in one go during the initialization phase and performs page resident, so that the communication between the central processing unit and the graphics processor is only through scalar signals of no more than 32 bytes during the iterative solution phase. This reduces the amount of data transferred per iteration on the peripheral component interconnect bus from tens of megabytes to tens of bytes of constant level, which is independent of the matrix size. This eliminates the constraint of the data transport bandwidth bottleneck on the solution throughput in the existing heterogeneous computing architecture.

[0018] Under the constraint matrix dimension of 50,000 rows by 200,000 columns, the latency of data movement in a single iteration is reduced from milliseconds to microseconds, and the utilization rate of the graphics processing unit (GPU) is increased from 20% to 30% to 85% to 95%. By employing a product-form basis inverse representation combined with a periodic reconstruction mechanism, explicit storage of the basis matrix inverse is avoided while maintaining numerical accuracy. This allows the method to handle large-scale sparse linear programming problems within the memory capacity limitations of a single GPU. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall process of a heterogeneous scientific computing acceleration method based on a zero-copy mechanism of video memory, provided in an embodiment of this application.

[0020] Figure 2 This is a schematic diagram of the data flow during the system initialization and video memory resident loading phases provided in an embodiment of this application.

[0021] Figure 3 This is a schematic diagram of scalar signal interaction between the central processing unit and the graphics processing unit during the iterative solution stage provided in an embodiment of this application.

[0022] Figure 4 The flowchart shows the forward solution process of the inverse basis representation of the product form provided in the embodiments of this application.

[0023] Figure 5 The flowchart shows the logic of the reverse solution process for the product form basis inverse representation provided in the embodiments of this application.

[0024] Figure 6 A schematic diagram of the tree-like reduction structure within a thread block for the two-level parallel reduction extremum algorithm provided in this application embodiment.

[0025] Figure 7 This is a schematic diagram illustrating the execution process of parallel binary search row by row in sparse compressed sparse row format column extraction provided in the embodiments of this application. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with the accompanying drawings and specific embodiments. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0027] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0028] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description of this application is provided in conjunction with the accompanying drawings and embodiments.

[0029] Method Implementation Examples This embodiment provides a heterogeneous scientific computing acceleration method based on a zero-copy memory mechanism. In a specific implementation, this method transmits the full constraint matrix of a large-scale linear programming problem to the high-bandwidth memory of the graphics processor (GPU) in a compressed sparse row format during the initialization phase and performs page-locked residency operations. During the iterative solution phase, the CPU performs the selection decision of entering and leaving base variables based solely on the scalar summary information returned by the GPU (no more than 32 bytes). The GPU autonomously completes all the intensive operations of parsimony-cost parallel computation, minimum ratio parallel verification, and implicit update of the base matrix within the high-bandwidth memory. This reduces the data transfer volume on the peripheral component interconnect bus for each iteration from tens of megabytes to a constant level of tens of bytes. This method solves the technical problem that frequent transfer of constraint matrix data between the CPU main memory and the GPU video memory in existing heterogeneous computing architectures causes bus bandwidth to become a bottleneck in solution throughput. It achieves the beneficial effect of increasing the GPU computing unit utilization from 20% to 30% to 85% to 95% under the condition that the constraint matrix has a dimension of 50,000 rows by 200,000 columns, and increasing the end-to-end solution throughput by four to six times.

[0030] In one specific embodiment, this method operates on a heterogeneous computing all-in-one machine equipped with an x86 architecture CPU and an NVIDIA architecture GPU. The CPU and GPU are connected via a 4th generation peripheral component interconnect bus (PCIe) with a 16-channel configuration. The theoretical bidirectional bandwidth of this bus is 32 gigabytes per second (32 GB / s), and the measured effective bandwidth is 12-15 gigabytes per second (12-15 GB / s). The GPU is equipped with high-bandwidth Gen 2 enhanced video memory (VRAM) with an on-chip access bandwidth of 2039 gigabytes per second (2039 GB / s) and a total capacity of 80 gigabytes (80 GB). The linear programming problem handled by this method is in standard form: solving for the objective function... Decision variable vector that achieves minimum value The constraints are and ,in For the constraint matrix, The right-hand constraint vector, This represents the cost vector for the objective function. In the large-scale sparse scenario addressed in this embodiment, the number of constraint rows... Not less than Number of decision variables Not less than constraint matrix The proportion of non-zero elements is less than one percent.

[0031] like Figure 1 As shown, this method includes the following steps: S000: The constraint matrix, right-hand constraint vector, and objective function cost vector of the linear programming problem are transferred from the CPU main memory to the high-bandwidth video memory of the graphics processor in one go. A page-locking operation is performed on the video memory region in the high-bandwidth video memory that stores the constraint matrix, so that the constraint matrix resides in the high-bandwidth video memory during subsequent iterative solutions.

[0032] like Figure 2 As shown, the data flow in this stage sequentially goes through three stages: format conversion and page-locked memory allocation at the central processing unit, one-time direct memory access transmission via the peripheral component interconnect bus, and video memory writing and page-locked residency at the graphics processing unit.

[0033] S010: Data import and format verification for linear programming problems.

[0034] The data acquisition interface on the central processing unit receives the raw data of the linear programming problem provided by the application, the raw data including the constraint matrix. Right-hand constraint vector and the objective function cost vector The constraint matrix The data is stored in the CPU main memory as a list of triples in coordinate format or a dense row-major array. Data integrity verification includes checking the constraint matrix. Is the number of rows equal to the right-hand constraint vector? Dimensions Check the constraint matrix Is the number of columns equal to the objective function cost vector? Dimensions And check the right-hand constraint vector Are all component values ​​not less than zero? When the constraint matrix... When inputting in coordinate format, a conversion operation is performed from coordinate format to compressed sparse row format—all triples are sorted by row index, generating an array of row pointers, an array of column indices, and an array of values.

[0035] The length of the row pointer array is (Each element is a 32-bit integer), the row pointer array of the... The element records the first... The starting position of the first non-zero element of the row in the value array. The length of the column index array is... (Each element is a 32-bit integer), where The constraint matrix The total number of non-zero elements in the array. The length of the value array is also... (Each element is a 64-bit double-precision floating-point number). The output of this step is a constraint matrix of three arrays and a dense vector stored in the CPU's main memory in a compressed sparse row format. and And the problem size parameter triplet .

[0036] For example, suppose the constraint matrix dimension of a small-scale verification instance is . , constraint matrix For one Given a matrix whose first row, first column contains 1, first row, third column contains 1, second row, second column contains 2, second row, fourth column contains 1, third row, first column contains 3, third row, second column contains 2, third row, fifth column contains 1, fourth row, first column contains 1, fourth row, second column contains 1, fourth row, sixth column contains 1, and all other elements are zero. What is the total number of non-zero elements in this matrix? After format conversion, the row pointer array is... The column index array is The value array is All elements are dimensionless values. Right-hand constraint vector. Cost vector In the row pointer array, the 0th element with a value of 0 indicates that the non-zero elements in the first row start from the 0th position of the value array; the 1st element with a value of 2 indicates that the non-zero elements in the second row start from the 2nd position, and so on. Verifying the third row: the 2nd element of the row pointer array is 4 and the 3rd element is 7, indicating that the non-zero elements in the third row occupy positions 4 to 6 in the value array, corresponding to column indices of... Value ,Right now , , It is consistent with the original matrix.

[0037] S020: Graphics processor device detection and video memory capacity pre-check.

[0038] The central processing unit calls the graphics processing unit device attribute query interface to obtain the total high-bandwidth video memory capacity of the target graphics processing unit. (Unit: bytes). Based on the aforementioned problem size parameter. The total amount of video memory required by this method is estimated according to the following. (Unit: bytes): The value array occupies bytes, the row pointer array occupies bytes, the column index array occupies Byte, right-hand constraint vector Occupy Bytes, cost vector Occupy Bytes, current basic feasible solution vector Occupy Bytes, simplified cost vector Occupy Byte, workspace buffer Let it be 10% of the sum of the aforementioned terms. In this embodiment, the term introduced by the Big M method... Each artificial variable is a standard basis vector with a fully defined structural characteristic, and does not need to be explicitly stored in a compressed sparse row format. When Exceed When 85% is reached, an insufficient memory error is reported and execution terminates. When the threshold is not exceeded, a memory allocation plan table is generated. This plan table records the starting offset address of each data block in the high-bandwidth memory of the graphics processor, including the offset address of the value array, the offset address of the row pointer array, the offset address of the column index array, the offset address of the right-hand constraint vector, the offset address of the cost vector, the offset address of the basic feasible solution, the offset address of the simplified cost vector, and the offset address of the workspace. All offset addresses are aligned to 8-byte boundaries to meet the natural alignment requirements of 64-bit double-precision floating-point types, avoiding performance loss caused by unaligned access.

[0039] For example, for the constraint matrix dimension , , Large-scale instances, the value array occupies The row pointer array occupies approximately 38.15 megabytes (MB). The column index array occupies approximately 0.19 megabytes (MB). The right-hand constraint vector occupies approximately 19.07 megabytes (MB). The cost vector occupies approximately 0.38 megabytes (MB). The basic feasible solution vector occupies approximately 1.53 megabytes (MB). Bytes, approximately 0.38 megabytes (MB), simplified cost vector occupancy. Bytes, approximately 1.53 megabytes (MB), the sum of the aforementioned items is This is approximately 61.24 megabytes (MB), and the workspace is taken as 10% of this value. That is approximately 6.12 megabytes (MB). Total approximately That's approximately 67.36 megabytes (MB), far smaller than the NVIDIA A100 graphics processor. 85% of gigabyte (GB) high-bandwidth video memory capacity Gigabytes (GB) meet the capacity pre-check requirements.

[0040] S030: Central Processing Unit Page Lock Memory Allocation and Data Population.

[0041] The central processing unit (CPU) calls the page-locked memory allocation interface of the graphics processing unit (GPU) runtime library to allocate memory on the CPU side in a size equal to... The page-locked memory region. Physical memory pages within this region are marked as non-swapable by the operating system, preventing subsequent direct memory access transfers from triggering virtual memory page swapping and ensuring that the transmission bandwidth consistently reaches the peak level of the peripheral component interconnect bus link. The compressed sparse row format array of row pointers, column indexes, and values, along with the right-hand constraint vector, are used to... and cost vector The data is copied sequentially to a contiguous address space within the page-locked memory region. The output of this step is a contiguous block of data already filled in the page-locked memory of the central processing unit.

[0042] S040: One-time direct memory access transfer to the graphics processor's high-bandwidth video memory and page-locked resident.

[0043] Calling the synchronous memory copy interface of the graphics processor runtime library, the entire page-locked memory region on the CPU side is locked. Byte data is transferred in a host-to-device manner to the base address of the high-bandwidth video memory of the graphics processor (GPU). This one-time direct memory access transfer refers to a batch transfer operation completed through a single synchronous memory copy interface call. The underlying hardware may break it down into multiple transaction packets according to the peripheral component interconnect bus protocol, but it appears as a single atomic operation to the caller. This transfer operation is executed through the direct memory access engine of the peripheral component interconnect bus. After the transfer is complete, the video memory access suggestion interface is called to mark the region in the high-bandwidth video memory containing the constraint matrix data as read-only priority access mode (i.e., setting this region to read-primarily access mode through the video memory suggestion interface), to guide the GPU's cache controller to adopt a more aggressive cache prefetch strategy for this region. Subsequently, the asynchronous video memory prefetch interface is called to prefetch the data pages in this region into the GPU's L2 cache to reduce cache miss latency during the first access in subsequent iterations. This transfer is the only large-scale data movement from the CPU main memory to the GPU video memory in this method.

[0044] For example, for the aforementioned An instance of approximately 67.36 megabytes (MB) achieved a measured effective bandwidth of approximately 13 gigabytes per second (13 GB / s) in a 16-channel configuration of the fourth-generation peripheral component interconnect bus standard. The time taken for this single transfer was approximately... This is approximately 5.2 milliseconds (ms). This transmission time is a one-time fixed overhead and will not occur again in subsequent thousands to tens of thousands of iterations.

[0045] S050: Initial basic feasible solution and metadata initialization.

[0046] This method employs the Big M method to construct the initial basic feasible solution. A parallel computing kernel is launched on the graphics processing unit to perform computation for each constraint row. ( from arrive Logically introduce an artificial variable The cost coefficient of the artificial variable is set as follows: The selection criterion is that this value must be much larger than the cost vector. The factor is 100 times the absolute value of any component in the equation. This factor ensures that the cost of the artificial variable is always positive within the normal iteration range, thus gradually driving it out of the basis during the optimal solution process. This is because the constraint matrix columns corresponding to the artificial variable... For the first standard basis vectors (No. Behavior 1, all others 0), possesses fully deterministic structural characteristics, requiring no explicit storage in the compressed sparse row format. All computational kernels on the graphics processor side that need to access matrix columns have column indices greater than […]. The standard basis vector is returned directly via a conditional branch. Initial basis variable index set. Stored in a 32-bit integer array at the central processing unit, with a length of The initial value is Non-basic variable index set Stored in another 32-bit integer array at the central processing unit, with a length of The initial value is Initial basic feasible solution It is already in place in the high-bandwidth video memory of the graphics processor. The current objective function value is initialized at the CPU level. Simultaneously, a scalar signal buffer of 32 bytes, containing four 64-bit fields, is allocated on the graphics processor side to store the summary data returned by the graphics processor to the central processing unit (CPU) in each iteration. A page-locked scalar receive buffer of the same size is also allocated on the CPU side.

[0047] For example, for the aforementioned , Small-scale validation instance, initial basic variable index set (Column indexes of artificial variables, dimensionless), set of indexes for non-basic variables (Column indices of original and slack variables, dimensionless). To simplify the verification derivation, the following small-scale example uses slack variables to directly construct the initial basis to avoid introducing artificial variables. This simplification does not affect the applicability of the general process of this method. Because in this example... For slack variables, It has been established, and slack variables can be used directly as the initial basis. Skipping the Big M method, initial basic feasible solution (Each component is dimensionless), initial objective function value The elementary transformation vector chain represented by the product form basis inverse is initially empty, and the chain length is... The system has entered the iteration-ready state.

[0048] S100: On the graphics processor side, the reduction cost vector operation is performed on all non-basic variables in parallel based on the constraint matrix residing in the high-bandwidth video memory, and the global minimum reduction cost value and its corresponding entering-basic column index are obtained from the reduction cost vector through parallel reduction.

[0049] like Figure 3As shown, during the iterative solution phase, the central processing unit (CPU) and the graphics processing unit (GPU) communicate only through scalar signals of no more than 32 bytes. After the GPU autonomously completes all intensive calculations in its high-bandwidth video memory, it sends the scalar summary information back to the CPU.

[0050] S110: Simplified cost parallel computing on the graphics processor side.

[0051] On the graphics processing unit (GPU), the cost reduction operation is performed in parallel for all non-basic variables. For each non-basic variable in the linear programming problem... Its simplified cost Defined as the original cost coefficient With dual variable vector Same constraint matrix List The difference of their inner products, i.e. The dual variable vector Obtaining this requires performing operations on the inverse of the basis matrix and the vector of cost coefficients of the basic variables. Product operation .

[0052] This method uses a product-form basis inverse representation to avoid requiring the inverse of the basis matrix. Explicit storage and computation. The product form basis inverse representation will... It can be expressed as a product of a series of elementary matrices: ,in Let be the inverse matrix of the initial basis. For the case using the Big M method, the initial basis is the identity matrix. Therefore . For the first The elementary matrix generated by the principal component transformation is only related to the identity matrix in the first principal component transformation. The column is different; the column is replaced with the elementary transformation vector. The product form basis inverse is represented in the high-bandwidth video memory of the graphics processor as a set of elementary transformation vectors. and its corresponding set of principal position indices Stored in the form of , where This represents the number of principal transformations performed.

[0053] like Figure 5 As shown, the dual variable vector Obtained through a reverse solution process. Specifically, let The initial value is equal to the cost coefficient vector corresponding to the current basic variable. Then, in the reverse order of the generation order of the elementary transformation vectors (i.e., ... Perform the following inner product assignment operations sequentially: Calculate With the elementary transformation vector inner product ,Will In the The value at each component position is updated to the inner product value, i.e. , The remaining components remain unchanged. The inner product operation is implemented on the graphics processor by calling the double-precision vector inner product interface of the linear algebra library, with a computational complexity of O(n) for each call. Floating-point operations. All. After the inner product assignment operation is completed, the result is That is The total computational complexity of this reverse engineering process is... This is a floating-point operation. The inner product assignment operation is mathematically asymmetrical to the scattering rank-one update operation in the forward solution process: each step of the forward solution modifies all components of the input column vector, while each step of the reverse solution modifies only a single component of the input row vector at the pivot position. This difference stems from the algebraic structure difference between left multiplication of an elementary matrix by a column vector and right multiplication of a row vector by an elementary matrix.

[0054] Obtain the dual variable vector Then, the graphics processor's parallel computing kernel is activated for each non-basic variable. Perform simplified cost calculation. The thread configuration of the parallel computing kernel is: grid dimension equal to There are 256 threads per block. Each thread is responsible for a non-basic variable. For simplified cost calculation, two cases need to be distinguished. When the non-basic variables... The column index is not greater than At that time, that is, the variable corresponds to the original decision variable column, and the thread traverses the constraint matrix. In the compressed sparse row format, the column index range for each row is located at the column index equal to... The non-zero elements and the dual variable vector Perform multiplication and summation operations on the corresponding row components, and the sum of the multiplication and summation results of all rows is the sparse inner product. Then from the cost vector Read from And perform subtraction to obtain When the non-basic variables The column index is greater than When, that is, the column of artificial variables corresponding to the variable, the column of the constraint matrix. Standard basis vectors Without needing to look up from the compressed sparse line format, the sparse inner product degenerates into a direct read. The Each component, cost coefficient Equal to the large M penalty coefficient The simplified cost is directly... The output of this step is a parsimonious cost vector stored in the high-bandwidth video memory of the graphics processor. .

[0055] For example, regarding the aforementioned , The small-scale verification instance executes step S110 of the first iteration. Initial base (Slack variables), the initial basis matrix is ​​the identity matrix. The product form of the basis inverse represents the chain length. Therefore, the reverse solution process does not perform any rank-one update. Directly equal to Non-basic variable index set .right : Column 1 of the constraint matrix , , .right : Column 2 of the constraint matrix , , Simplified cost vector .

[0056] S120: Parallel reduction on the graphics processor side to find the global minimum parsimony cost.

[0057] like Figure 6 As shown, from the simplified cost vector Obtain the global minimum parsimony cost value and its corresponding base column index It employs a two-stage parallel reduction structure to execute in parallel on the graphics processing unit. In the first-stage reduction, Each simplified cost value is evenly distributed to Each thread block contains 256 threads, each responsible for loading its own simplified cost value into the shared memory of the thread block. Tree-structured comparison reduction is performed in this shared memory: in each iteration, the number of active threads is halved, and each active thread compares the value at its current position with the value at the step position, retaining the smaller value and its global index. After each round of iterations, thread 0 of each thread block holds the local minimum value and its global index for that thread block, and writes the result to a local result array in global memory. In the second-level reduction, a single thread block is started to perform the same tree-structured comparison reduction on all local results to obtain the global minimum parsimony cost value. and the corresponding base column index The total number of comparisons for this two-level reduction is The parallel depth is... Level. The thread block size is fixed at 256, which is 8 times the NVIDIA graphics processor thread bundle size of 32. Under this configuration, the 32 shared memory banks do not experience access conflicts. Write the first 16 bytes of the scalar signal buffer on the graphics processor side.

[0058] For example, continuing from the first iteration mentioned above, The reduction can be completed with only one thread block. Thread 0 holds the... Thread 1 holds In the first round of tree-based reduction, thread 0 compares... and , ,reserve and its index . .

[0059] S130: Scalar return and optimality determination at the central processing unit.

[0060] The peripheral component interconnect bus scalar signal channel manager calls the asynchronous memory copy interface to copy 16 bytes of data (including...) from the graphics processor's scalar signal buffer. and The data is asynchronously copied from device memory to the CPU's page-locked scalar receive buffer, and then the stream synchronization interface is called to wait for the transmission to complete. The CPU reads the data... The value is compared with the negative numerical precision tolerance. Comparison, among which The tolerance is selected based on the machine precision of 64-bit double-precision floating-point. Take it This is to accommodate the floating-point rounding error accumulated during the iteration process. When When the current solution satisfies the optimality condition, and the reduction costs of all non-basic variables are non-negative, proceed to S400 to execute the convergence termination procedure. At that time, confirm For the column index of the base variable, continue executing S200.

[0061] For example, continuing from the first iteration mentioned above, , If true, then the current solution does not satisfy the optimality condition. Index for the base variable column.

[0062] S200: When the current solution does not satisfy the optimality condition, the constraint matrix column vector corresponding to the base column index is extracted from the high-bandwidth video memory at the graphics processor end, and the basis inverse product operation is performed on the column vector to obtain the updated column vector. Based on the updated column vector and the current basic feasible solution, the minimum ratio test is performed in parallel to obtain the off-base row index and the minimum ratio, and the off-base row index and the minimum ratio are sent back to the central processing unit in the form of a scalar signal.

[0063] S210: Graphics processor-side input basis column vector extraction and basis inverse product.

[0064] The in-base index is extracted from the high-bandwidth video memory at the graphics processor side. The corresponding constraint matrix column vectors need to be distinguished into two cases. For example... Figure 7 As shown, when The value is not greater than At that time, that is, when the entering variable corresponds to the original decision variable column, the parallel column extraction calculation kernel is started for each row of the constraint matrix. ( from arrive Allocate a thread for parallel execution. Each of the aforementioned threads is in the first... Compressed sparse row format column index range Perform a binary search to locate the column index equal to The element. The binary search process is as follows: initialize the lower bound of the search. and the upper bound of the search ;when When, take the middle position. ,like equal Then The first one written into the output vector Each component is returned, if... Less than Then update Otherwise update When the search interval is exhausted and no result is found. When, the first of the output vectors will be... Setting the nth component to zero indicates that the constraint matrix is ​​in the nth... Line 1 The columns contain zero elements. Because the compressed sparse row format guarantees that the column indices within each row are sorted in ascending order, the correctness of the binary search is ensured. The search complexity for each row is O(n log n). ,in For the first The number of non-zero elements in a row. The lines are executed simultaneously in parallel, with a total parallel depth of [value missing]. The thread configuration of the parallel column extraction computing kernel is as follows: Each thread block contains 256 threads. When... The value is greater than At that time, the column of artificial variables corresponding to the basic variables, and the column of the constraint matrix. For the first standard basis vectors The standard basis vectors are generated directly on the graphics processor side. Each component is The remaining components are There is no need to look up from the compressed sparse line format.

[0065] After the column vector extraction is completed, the extracted column vectors are... Perform the basis inverse product operation to obtain the updated column vector. .like Figure 4 As shown, this operation is achieved through a forward solution process. Let... The initial value is equal to Then, according to the order in which the elementary transformation vectors were generated (i.e. Perform the following rank-one update operations sequentially: Read In the The current value at each component location is recorded as a temporary scalar. ,Will The Clear all components to zero, then perform vector addition. The vector addition is implemented on the graphics processor by calling the double-precision vector scalar multiplication and addition interface of the linear algebra library. After the update of Step Rank 1 is completed, That is, the updated column vector The forward solution process differs from the reverse solution process in S110 in two ways: firstly, the iteration directions are different; the forward solution proceeds according to... to Execute in forward order, solve in reverse order to The first step involves reversing the order of operations. Secondly, the single-step operations differ. Each step in the forward solution is a scattering rank-one update operation—extracting the component value of the input column vector at the pivot position, clearing it to zero, and then adding the product of this value and the elementary transformation vector to all components. Each step in the reverse solution is an inner product assignment operation—calculating the inner product of the input row vector and the elementary transformation vector, and assigning the inner product value to the component at the pivot position while leaving the other components unchanged. The forward solution is used to obtain the product of the inverse of the basis matrix and the column vector. The reverse solution is used to obtain the product of the row vector and the inverse of the basis matrix. The total computational complexity of this forward solution process is... Floating-point operations.

[0066] For example, continuing from the first iteration mentioned above, , This corresponds to the original variable column. For the constraint matrix... Perform a binary search on each row to locate the column index. The elements (Note: internal column indices start from 0, (This corresponds to the second column, i.e., index 1). Taking the third row as an example: , The search range is . , , .initial , , , ,equal (0-indexed), hit. Extracting column vectors (Each component is dimensionless). The product form basis inverse represents the chain length. Forward solving does not perform any rank-1 update. (Each component is dimensionless).

[0067] S220: Minimum Ratio Check for Parallelism on the Graphics Processor Side.

[0068] Based on the updated column vector Compared with the current basic feasible solution Perform minimum ratio tests in parallel to obtain off-base row indexes and minimum ratio Start the parallel computing kernel and configure the threads as follows: Each thread block contains 256 threads. ( from arrive Perform the following operations: Check Is it greater than the principal component numerical stability threshold? The threshold was chosen based on 64-bit double-precision floating-point machine precision. of This is to avoid numerical instability caused by dividing by values ​​close to zero. At that time, the thread calculates the ratio. .when Not greater than At that time, the thread will Let it be positive infinity.

[0069] The global minimum value and its index of the ratio vector are obtained through the same two-level parallel reduction structure as S120. In the degenerate case where the ratios of multiple basic variable rows are all equal to the minimum ratio, the Brand anti-loop rule is applied: in all... Select the basic variable index set in the row. The row with the smallest index value is used as the off-base row index. This rule theoretically guarantees that the simplex method terminates within a finite number of steps without falling into a degenerate cycle. Write the last 16 bytes of the scalar signal buffer on the graphics processor side.

[0070] For example, continuing from the first iteration mentioned above, (Each component is dimensionless) (Each component is dimensionless). For the first... OK: , For the first OK: , For the first OK: , For the first OK: , . , (Dimensionless, corresponding to the position of the second basic variable).

[0071] S230: Scalar return and off-base confirmation at the central processing unit.

[0072] The peripheral component interconnect bus scalar signal channel manager performs the same 16-byte asynchronous direct memory access passback operation as S130, transferring the data from the graphics processor's scalar signal buffer... Copy from device memory to the CPU's page-locked scalar receive buffer. The CPU reads the... The value of . When When equal to positive infinity (implementation check) If the linear programming problem is determined to be unbounded, an unbounded error signal is output and execution is terminated. When the value is finite, confirm For off-base row indexes, perform CPU-side metadata update operations: index arrays in non-basic variables. Search in The position, replace the value at that position with The original value, then indexed in the base variable array. Lieutenant General Each element is updated to Update the current objective function value. The output of this step is the updated set of basic variable indices on the central processing unit. Non-basic variable index set and objective function value and off-base row index .

[0073] For example, continuing from the first iteration mentioned above, , , . The original value is (slack variables) (column index), will In Replace with ,Will Updated to After the update , . .

[0074] S300: On the graphics processor side, perform an implicit update operation of the basis matrix in-situ in the high-bandwidth video memory based on the advancing column index and the departing row index, and update the current basic feasible solution.

[0075] S310: Receive master element transformation instruction.

[0076] The central processing unit will use the base column index. and the off-base row index The data is packaged into 8 bytes (two 32-bit integers) and sent from the CPU-side page-locked memory to the GPU-side instruction buffer via the peripheral component interconnect bus scalar signal channel manager, using the asynchronous memory copy interface. This 8-byte transfer occurs over the peripheral component interconnect bus, primarily due to direct memory access transaction establishment latency, and lasts approximately 1 to 2 microseconds.

[0077] S320: Implicit update operation of the basis matrix, namely, elementary transformation vector generation and pivot row normalization.

[0078] On the graphics processor side, based on the updated column vector and the off-base row index Generate the elementary transformation vector corresponding to this principal component transformation. The pivot row and column transformation described in this application refers to an implicit basis matrix update operation implemented by generating elementary transformation vectors and appending them to the transformation chain of the product form basis inverse representation, rather than performing an explicit Gaussian elimination row and column transformation on the constraint matrix. The pivot value of the updated column vector at the position furthest from the basis row index is read. Obtain the reciprocal of the pivot value. Start the parallel computing kernel and configure the threads as follows: Each thread block contains 256 threads. For each line... ( from arrive ),when When, the first of the elementary transformation vectors Each component is ;when When, the first of the elementary transformation vectors Each component is The elementary transformation vector The first elementary transformation vector array stored in the high-bandwidth video memory of the graphics processor in product form, represented by the basis inverse, is... Each position records the principal position index. Update chain length The elementary transformation vector array has a pre-allocated capacity of [amount] in the high-bandwidth video memory of the graphics processor. bytes, of which This is the preset reconstruction threshold.

[0079] The elementary transformation vectors described in this method The dimensional system is explained as follows: The constraint matrix The elements are coefficients of the linear programming problem, which are dimensionless numerical values; the updated column vector It is the product of the inverse of the basis matrix and the column vectors of the constraint matrix. Each element in the inverse of the basis matrix is ​​also a dimensionless value. Each component is a dimensionless value. (The reciprocal of the principal component...) The reciprocal of the dimensionless scalar is also a dimensionless numerical value. The elementary transformation vector... Each component is a dimensionless value, and the vector scalar multiplication and addition operations are performed in both the forward and reverse solving processes. middle, It is a dimensionless scalar. It is a dimensionless vector. The dimensions depend on the physical meaning of the input vector (such as the cost coefficient vector). A dimensionless, fundamentally feasible solution (The values ​​of the decision variables are also dimensionless).

[0080] For example, continuing from the first iteration mentioned above, (Each component is dimensionless) principal component value , The elementary transformation vector The components are: , (Principal position, dimensionless). , The elementary transformation vector (Each component is dimensionless), principal component position Chain length .

[0081] S330: Basic feasible solution update.

[0082] On the graphics processor side, based on the updated column vector The minimum ratio and the off-base row index Update the current basic feasible solution stored in high-bandwidth video memory. Start the parallel computing kernel and configure the threads as follows: Each thread block contains 256 threads. For each line... ( from arrive ),when When, execute ;when When, execute This update operation ensures that the new basic feasible solution satisfies the feasibility conditions. :for OK, exist The time obviously does not decrease, in The minimum ratio test guarantees that thus .

[0083] For example, continuing from the first iteration mentioned above, , , , . , (dimensionless, Assignment), , After the update .

[0084] S340: Basic variable index graphics processor self-updates.

[0085] On the graphics processor side, based on the principal transformation instructions received in step S310 (including the base column index) and off-base row index The graphics processor (GPU) maintains a copy of the locally maintained basic variable index set, which it updates automatically. The GPU also maintains a copy of the basic variable index array in high-bandwidth video memory. Local copy (size) (bytes), after the local copy is synchronized from the central processing unit to the graphics processing unit (GPU) via a one-time direct memory access transfer during the S050 initialization phase, it will not receive full updates of the array from the central processing unit during subsequent iterations. After each pivot transformation, the GPU starts a single-threaded computing kernel to perform an atomic index update operation: The value is updated to ,Right now This operation involves only one 32-bit integer memory write, completed locally within the high-bandwidth video memory of the graphics processor, without generating any data transfer on the external component interconnect bus. Therefore, during the iteration phase, the amount of data transferred between the CPU and the graphics processor per iteration is strictly limited to scalar digest return (no more than 32 bytes) plus principal component transformation instruction transmission (8 bytes), totaling no more than 40 bytes, regardless of the size of the constraint matrix. The CPU-side basic variable index array... Local copy with graphics processor By performing the same update operation Maintain logical consistency without explicit synchronization.

[0086] S350: Product form basis inverse representation of periodic reconstruction.

[0087] After each execution of S320, check the cumulative number of elementary transformation vectors in the product form basis inverse representation. Does it exceed the preset reconstruction threshold? The preset reconstruction threshold The range of values ​​is For example, it can be set to The selection of this threshold is based on the following technical teaching: each call to the forward and backward solving of the inverse representation of the product form's basis requires execution... Each of the vector scalar multiplication and addition operations has a computational complexity of O(n log n). Therefore, the total complexity of a single-radix inverse product is (number of floating-point operations). .when The computation time exceeds that of a single triangular decomposition of a basis matrix. When this occurs, a refactoring should be triggered to reset the chain length. For the NVIDIA A100 graphics processor (equipped with 108 streaming multiprocessors), in Under these conditions, a single vector scalar multiplication and addition operation performed using a linear algebra library takes approximately [time missing]. A single triangular decomposition takes approximately [time in microseconds] to execute via the solver library. milliseconds, theoretical upper limit However, considering the floating-point rounding error of elementary transformation vector chain product varies... Linear accumulation, in The time-product error reaches the order of magnitude This may lead to incorrect sign determination for simplified expenses, therefore the actual amount taken... As a conservative choice to balance computational efficiency and numerical accuracy. In some cases, the fixed overhead of frequent triangular decomposition actually outweighs the computational cost saved in chain product operations.

[0088] when At that time, on the graphics processor side, the current basis matrix is... Perform partial principal component triangulation. Specifically, call the double-precision triangulation interface of the graphics processor solver library. Doolittle decomposition of line swapping ,in For permutation matrix, It is a lower triangular matrix. It is an upper triangular matrix. The triangular decomposition is parallelized on the graphics processor through block matrix multiplication. After decomposition, based on... and Factors are regenerated into product form, basis inverse representations are added, and chain lengths are reset. .

[0089] This method also includes an adaptive accuracy monitoring mechanism. A preset verification interval is executed each time. After the next iteration, a lightweight verification calculation kernel is launched on the graphics processor: indexing the set of basic variables. Select a basic variable column index (Example, selection is made sequentially using a round-robin method or using a uniform pseudo-random method), and the forward solution process is used to... The corresponding constraint matrix column vectors are subjected to basis inverse product operation to obtain the verification vector. If the numerical accuracy of the basis inverse representation of the product form is good, then It should be the first standard basis vectors The infinite norm residual between the verification vector and the expected standard basis vector. Same as preset accuracy threshold When comparing, At that time, the verification calculation kernel writes the reconstruction trigger flag into a pre-allocated single-byte flag buffer in the high-bandwidth video memory of the graphics processor and terminates execution. The central processing unit (CPU) reads the flag buffer through synchronous memory copying after each verification interval. When the flag is set, the CPU calls the double-precision triangulation interface of the graphics processor solver library to perform the aforementioned partial principal component triangulation to restore numerical precision. After completion, the flag is cleared to zero. This adaptive precision monitoring mechanism trades approximately 5% additional computational overhead for numerical stability assurance.

[0090] The iterative cycle is as follows: S100, S200, and S300 are executed sequentially to form a complete simplex method iteration, and the cycle continues to be executed until the optimality condition is determined in S130.

[0091] S400: Convergence Termination and Result Feedback.

[0092] When S130 determines the global minimum simplified cost value Not less than the negative numerical precision tolerance When the convergence termination process is executed, proceed.

[0093] S410: Feasibility check of artificial variables.

[0094] The central processing unit traverses the array of base variable indices. All For each element, check if there exists an index value greater than 1. Artificial variables remain in the position of basic variables. When a basic variable exists... The value is greater than And its corresponding basic feasible solution components Greater than the preset feasibility tolerance When the original linear programming problem has no feasible solution, it is determined that the original linear programming problem has no feasible solution. When there exists an artificial variable in the basis but its corresponding basic feasible solution component value is not greater than... When this condition is met, it is considered a degenerate feasible solution and thus feasible. When no artificial variables exist in the basis, it is directly deemed feasible. The computational complexity of this traversal is O(n log n). It takes about tens of microseconds on the central processing unit.

[0095] S420: Optimal solution vector backhaul.

[0096] Construct the complete decision variable solution vector on the graphics processor. The parallel computing kernel is activated, allocating a memory slot in the high-bandwidth video memory of the graphics processor. A 3D double-precision floating-point array initialized to an all-zero vector. Each basic variable row ( from arrive The thread will Write the value The first vector Each location (scattering operation). Then, via a synchronous memory copy interface... Vector (size) (Bytes) are transferred from the graphics processing unit's device memory to the central processing unit's main memory in a single operation. This is the only large-scale data transfer from the graphics processing unit to the central processing unit in this method. The output is the optimal solution in the central processing unit's main memory. and optimal objective value .

[0097] S430: Release video memory resources.

[0098] The system calls the graphics processor runtime library's memory release interface to release all allocated data blocks in the graphics processor's high-bandwidth memory sequentially, including the constraint matrix array, the elementary transformation vector array in product form (basis inverse representation), the pivot index array, the scalar signal buffer, the instruction buffer, the basic variable index buffer, the working vector, and the temporarily allocated solution vector. The system then calls the page-locked memory release interface to release all page-locked memory regions on the CPU side. System resource reclamation is complete.

[0099] For example, continuing the complete solution process of the aforementioned small-scale verification example, after the first iteration, the results have been obtained. , , , , , .

[0100] Complete derivation of the second iteration: Execute S110, Non-basic Variable Index Set .right :extract (Each component is dimensionless), call forward solver. , . : , ,because Vectors remain unchanged. . (Each component is dimensionless). Current (Each component is dimensionless). .

[0101] right : (Each component is dimensionless), call forward solver. , . : , Clear have to ,implement . (Each component is dimensionless). .

[0102] Execute S120, , , .

[0103] Execute S130, The optimality condition is not met.

[0104] Execute S210, (Each component is dimensionless).

[0105] Execute S220, . , ( ), , . A draw would occur if the score is tied. Applying the Brand anti-loop rule, , , Select .

[0106] Execute S230, Original value ,Will In Replace with , Updated to . , . .

[0107] Execute S320, , principal component value , . , , (Principal position, dimensionless). . , , .

[0108] Execute S330, , , ( Assignment (dimensionless). . .

[0109] Complete derivation of the third iteration: Execute S110, currently , , , .

[0110] right : Slack variable column (Each component is dimensionless). Call the forward solver. , . : , , constant. : , Clear have to , (Each component is dimensionless). . .

[0111] right : Slack variable column (Each component is dimensionless). Call the forward solver. , . : , Clear have to , (Each component is dimensionless). : , Clear have to , (Each component is dimensionless). . .

[0112] Execute S120, , , .

[0113] Execute S130, The optimality condition is met. Jump to S400.

[0114] Execute S410, traverse No artificial variables (all indexes) ) in the base.

[0115] Execute S420, ,according to Mapping: , , , Non-basic variables , . (Each component is dimensionless) .

[0116] Peripheral component interconnect bus transmission statistics: One-time transmission during initialization phase The data volume is 32 bytes (including constraint matrix data and an initial copy of the basic variable index array). During each iteration, a 32-byte scalar digest is sent back, plus an 8-byte principal component transformation instruction. The graphics processor autonomously updates its local basic variable index set copy based on these instructions without incurring additional transmission. The bus data transfer volume for two iterations is... Bytes. Termination phase feedback. The total data transfer volume is approximately 348 bytes (excluding the one-time transfer during the initialization phase). Traditional methods require transferring the entire base matrix column data in each iteration, with the transfer volume being the dimension of the matrix column multiplied by 8 bytes multiplied by the required number of columns. In large-scale problems, this difference can be magnified to several orders of magnitude.

[0117] This embodiment provides a system for executing the aforementioned heterogeneous scientific computing acceleration method based on a zero-copy memory mechanism. The system runs on a heterogeneous computing all-in-one machine equipped with a central processing unit (CPU) and a graphics processing unit (GPU), which are connected via a peripheral interconnect bus (PEBBus). The system consists of five functional modules: a matrix preprocessing and memory-resident loading engine, a lightweight instruction decision-maker on the CPU side, a fully resident parallel principal component computing core on the GPU side, a PEBBus scalar signal channel manager, and a convergence determination and result return terminator. In practice, these five functional modules can be implemented by the CPU and GPU respectively executing instructions stored in their respective memories.

[0118] The matrix preprocessing and memory-resident loading engine runs on the main thread of the central processing unit and is responsible for performing one-time data preparation and transfer during system startup. This engine comprises three subunits: a format conversion subunit, a memory planning subunit, and a direct memory access transfer subunit. The format conversion subunit is configured to perform the conversion from coordinate format to compressed sparse row format, achieved through radix sorting and row pointer prefix sums, with a sorting complexity of O(log n). The memory planning subunit is configured to estimate the offset address of each data block in the high-bandwidth memory of the graphics processor, ensuring 8-byte alignment to meet the natural alignment requirements of 64-bit double-precision floating-point types. The direct memory access transfer subunit is configured to call the graphics processor runtime library interface to perform large-block transfers from the host to the device, utilizing all channels of the peripheral component interconnect bus to achieve an effective bandwidth of approximately 12 to 15 gigabytes per second. The core data structures operated by this engine include: a memory allocation schedule table, a structure containing eight 64-bit unsigned integer fields (72 bytes in total), recording the starting offset address of each data block, stored in the CPU stack memory; and a compressed sparse row format three-array, including a row pointer array ( (32-bit integer elements), column index array ( (32-bit integer elements) and value arrays ( Each 64-bit floating-point element is stored in the CPU's page-locked memory and then transferred to the GPU's high-bandwidth video memory. The engine exits after initialization and has no runtime communication with the fully resident parallel principal component computing core.

[0119] The lightweight instruction decision unit at the central processing unit (CPU) runs on a single core of the CPU, performing all operations as scalar operations and small array indexing operations, and occupies approximately [memory value missing]. Bytes. For , For example, this decision-maker uses approximately [amount missing] memory. Megabytes (MB). The decision-maker comprises three subunits: an optimality determination subunit, an unboundedness determination subunit, and a basic variable index update subunit. The optimality determination subunit is configured to return the global minimum parsimony cost value and tolerance from the graphics processor. Perform a single floating-point comparison operation. The unboundedness determination subunit is configured to compare the minimum ratio with a positive infinity threshold. A comparison is performed. The basic variable index update subunit is configured to perform element swaps between the basic variable index array and the non-basic variable index array—a linear lookup of the in-basic column index in the non-basic variable index array. Location (complexity) (This process takes approximately tens of microseconds on a central processing unit, far less than the computation time on a graphics processing unit). Replace the value at that position with the original column index of the off-basic variable, and then... Each element is updated to The core data structure for this decision-maker operation includes: a basic variable index array. ( (32-bit integer elements), non-basic variable index array ( (32-bit integer elements), a page-locked scalar receive buffer (32 bytes, 4 64-bit floating-point elements), and the current objective function value. (64-bit floating-point scalar). The decision-maker receives the scalar digest returned by the graphics processor through the scalar signal channel manager, and sends an 8-byte principal transformation instruction (including the in-basic column index and the out-basic row index) to the graphics processor through the same channel.

[0120] The fully resident parallel principal component computing core is deployed on all the streaming multiprocessors of the graphics processor. Taking the NVIDIA A100 graphics processor as an example, this model is equipped with 108 streaming multiprocessors, each containing 64 parallel computing cores, for a total of 6912 parallel computing cores; the high-bandwidth video memory capacity is 80 gigabytes, and the on-chip bandwidth is 2039 gigabytes per second; each streaming multiprocessor is equipped with 48 kilobytes of shared memory. All computing cores read and write data from the high-bandwidth video memory without the involvement of the peripheral component interconnect bus. The computing core contains five sub-cores: a simplified cost parallel computing core, a minimum ratio parallel verification core, a basis matrix implicit update core (i.e., a principal component row and column transformation core), a product form basis inverse forward solver, and a product form basis inverse backward solver. The simplified cost parallel computing core is configured to process each non-basic variable... Perform sparse inner product operations with minimal cost, including manual variable decision branches. Thread configuration assigns one thread to each non-basic variable, with a thread block size of 256. The minimum ratio parallel check kernel configuration is used for... The ratio is calculated in parallel and two-level reduction is performed. The principal component row and column transformation kernel is configured to generate and store elementary transformation vectors and update the basic feasible solution. The forward solver is configured to perform calculations in parallel. to The linear algebra library interface is called in forward order to perform vector scalar multiplication and addition operations, used for updating column vectors in S210. The inverse solver is configured to... to The same operation is performed in reverse order and is used for the calculation of the dual variable vector in S110.

[0121] The core data structures for this computational core operation all reside in the high-bandwidth video memory of the graphics processor, including: a three-array constrained matrix compressed sparse row format (read-only, never modified from initialization to termination), and an array of elementary transformation vectors in product form with basis inverse representation (pre-allocated). Byte, for , Approximately 40 megabytes) and the principal position index array ( (32-bit integer elements, 400 bytes) The working vector includes the basic feasible solution. ( (64-bit floating-point elements), simplified cost vector ( (64-bit floating-point elements), dual variable vector ( (64-bit floating-point elements) and update column vector ( (64-bit floating-point element).

[0122] The Peripheral Component Interconnect Bus Scalar Signal Channel Manager is a thin logical layer running within the graphics processor runtime library on the CPU side, responsible for managing asynchronous direct memory access (DMI) transfers between the CPU and the graphics processor. This channel manager contains two sub-channels: an instruction downlink channel (configured for transmitting 8-byte major transformation instructions, host-to-device direction) and a digest uplink channel (configured for transmitting 32-byte scalar digest data, device-to-host direction). Since each transfer is only 8 to 32 bytes, actual throughput is dominated by the DMI transaction establishment latency (approximately 1 to 2 microseconds per transfer), rather than bus bandwidth limitations. The graphics processor autonomously performs atomic update operations on its local base variable index set copy based on the entering base column index and leaving base row index carried in the major transformation instruction, without requiring full synchronization of the base variable index array via the bus. This channel manager uses the graphics processor streaming mechanism to implement asynchronous non-blocking transfers, ensuring data consistency through a streaming synchronization interface. The instruction buffer is stored in two 8-byte copies (one in CPU-side page-locked memory and one in graphics processor-side device memory), containing two 32-bit integer fields: the entering base column index and the leaving base row index. The scalar signal buffer is stored in 32-byte duplicates and contains four 64-bit fields: global minimum parsimony cost value, in-basic column index, minimum ratio, and out-basic row index.

[0123] The convergence determination and result feedback terminator are logic modules. The convergence determination is performed in the instruction decision unit at the central processing unit (scalar comparison operation), and the result feedback is performed once through the scalar signal channel manager. This is a direct memory access transfer from the device to the host side, involving bytes. The finalizer consists of three subunits: a manual variable checking subunit (which iterates through the base variable index array on the CPU side, with a complexity of...). The solution vector assembly kernel (a parallel scattering operation on the graphics processor side, mapping sparse basic feasible solutions to basic variable indices) (The complete solution vector) and the return sub-unit (calling the synchronous memory copy interface to return the solution vector).

[0124] It should be noted that, in addition to using the aforementioned product form of basis inverse representation, the basis inverse product operation can optionally be implemented using the Bartels-Golub update method. The Bartels-Golub update method maintains the triangular decomposition of the basis matrix. Each time a principal component transformation occurs, it is updated via Givens rotation or Householder transformation. Factors and Factoring avoids complete refactoring. This alternative method has better numerical stability than the product form's basis inverse representation because it avoids the accumulation of floating-point rounding errors in long-chain products. However, the computational complexity of each update is O(n). Floating-point operations (Givens rotation needs updating) (A row of factors), higher than the product form basis inverse representation Parallel Givens rotations can be implemented on a graphics processor by calling the orthogonalization interface of the solver library. This alternative method is suitable for scenarios requiring high numerical accuracy but with a limited number of iterations.

[0125] In addition to employing the aforementioned row-by-row parallel binary search method, the column vector extraction can optionally include pre-constructing redundant copies of the compressed sparse column format during the initialization phase. These redundant copies of the compressed sparse column format enable column extraction from... The complexity of binary search is reduced to The direct access complexity, where For the first The number of non-zero elements in the column, but requires additional... Bytes of video memory. For For example, it provides approximately 600 megabytes (MB) of additional video memory, making it a better choice in scenarios with ample video memory.

[0126] In addition to employing the aforementioned custom computation kernel, the two-level parallel reduction can optionally directly call the device-level parameter minimum reduction interface of the parallel primitive library. Internally, this library uses a similar hierarchical reduction strategy but has undergone performance optimization for different graphics processor architectures, utilizing thread-beam-level data shuffling instructions to replace shared memory reduction. The advantage of using this library is its automatic adaptation to different graphics processor architectures; the disadvantage is the introduction of additional temporary memory allocation overhead.

[0127] It should be noted that when the system detects continuous The objective function value in the next iteration Change All less than At this time, the central processing unit's instruction decision-maker outputs a degradation warning log but continues execution. The Brand anti-loop rule theoretically guarantees that the simplex method can escape the degradation loop within a finite number of steps in the degradation case. When continuous... If the objective function value still does not improve in the next iteration, it is determined that there is a numerical problem in the basis inverse representation of the current product form, and the periodic reconstruction of S350 is triggered to restore the numerical accuracy through partial principal component triangular decomposition.

[0128] For example, a degenerate linear programming problem is constructed to verify the anti-loop mechanism. Let... (Each element is dimensionless) (Each component is dimensionless) (Each component is dimensionless) , Initial base , (Each component is dimensionless). During the iteration process, if the maximum coefficient selection rule is used instead of the Brand rule, multiple consecutive occurrences of degenerate vertices may occur at the degenerate vertex. If the principal component transformation does not improve the objective function value, the iteration will not terminate. Step S220 of this method enforces the Brand anti-loop rule, ensuring termination within a finite number of steps in all degenerate cases.

[0129] When the size of a linear programming problem exceeds the high bandwidth memory capacity of a single graphics processor (such as...) , , The required video memory is approximately 12 gigabytes, used only for compressing sparse row data; the total requirement, including the product form basis inverse representation and workspace, may exceed 80 gigabytes. Optionally, a block-based residency strategy is introduced. The constraint matrix is ​​partitioned column-wise into multiple blocks, each no larger than 70% of the available high-bandwidth video memory capacity of the graphics processor. The current iteration only needs to access the basis matrix columns and the entering basis columns, residing the partition containing these columns in video memory. When a change in the set of basis variables causes the basis column set to cross partitions, a small number of peripheral component interconnect bus transfers are triggered to swap in the new partition. Under this strategy, the data movement elimination rate is reduced from over 95% to approximately 70% to 85%, but it is still far superior to traditional solutions. In a multi-graphics processor configuration, different partitions can be assigned to different graphics processors to further increase the capacity limit.

[0130] The above is a detailed description of an implementation method for accelerating heterogeneous scientific computing based on a zero-copy memory mechanism provided in this application. It should be noted that the above embodiments are only used to illustrate the technical solution of this application and are not intended to limit the scope of protection of this application. For those skilled in the art, various modifications and changes can be made based on the above embodiments without departing from the spirit and scope of this application, and all such modifications and changes fall within the scope of protection of this application.

[0131] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit described above can be implemented in hardware.

Claims

1. A method for accelerating heterogeneous scientific computing based on a zero-copy memory mechanism, characterized in that, include: The constraint matrix, right-hand constraint vector, and objective function cost vector of the linear programming problem are transferred from the main memory of the central processing unit to the high-bandwidth video memory of the graphics processing unit in one go. The page-locking operation is performed on the video memory area storing the constraint matrix to keep the constraint matrix resident during the iterative solution process. The current basic feasible solution is initialized in the high-bandwidth video memory based on the right-hand constraint vector. On the graphics processor, a reduction cost operation is performed on all non-basic variables based on the resident constraint matrix to obtain the global minimum reduction cost value and the corresponding entering basis column index. This value is then sent back to the central processing unit in the form of a scalar signal to determine whether the current solution satisfies the optimality condition. When the optimality condition is not met, the column vector corresponding to the base column index is extracted at the graphics processor and the basis inverse product operation is performed to obtain the updated column vector. Based on the updated column vector and the current basic feasible solution, the minimum ratio test is performed to obtain the off-base row index. The off-base row index is returned to the central processing unit in the form of a scalar signal. Based on the advancing column index and the departing row index, the product form of the basis inverse is used to perform implicit updates of the basis matrix in-situ within high-bandwidth video memory and update the basic feasible solution. The central processing unit updates the basic variable index set and the non-basic variable index set and sends instructions to the graphics processor to synchronize the copy of the basic variable index set. Then, it returns to perform parsimony cost operations until the optimality condition is met.

2. The method according to claim 1, characterized in that, The one-time transmission uses a compressed sparse line format, including: Perform a format conversion on the constraint matrix to generate a row pointer array, column index array, and value array in a compressed sparse row format; The total capacity of the high-bandwidth video memory is obtained by calling the graphics processor device detection interface, and the required video memory is estimated based on the total number of non-zero elements, the number of constraint rows, and the number of decision variables in the constraint matrix. When the required video memory exceeds a preset percentage threshold of the total capacity, the execution is terminated. The central processing unit allocates a page-locked memory region and fills the page-locked memory region with the row pointer array, the column index array, the value array, the right-hand constraint vector, and the objective function cost vector. All data in the page-locked memory region is written to the pre-planned offset address of the high-bandwidth video memory via a one-time direct memory access transfer.

3. The method according to claim 1, characterized in that, The implicit update operation of the base matrix includes: Obtain the pivot value of the updated column vector at the index position of the off-base row, and obtain the reciprocal of the pivot value; Based on the reciprocal, a linear transformation is performed on each row component of the updated column vector to generate an elementary transformation vector, and the elementary transformation vector is appended to the transformation chain of the basis inverse representation.

4. The method according to claim 3, characterized in that, The basis inverse product operation uses the basis inverse representation of the product form, which expresses the inverse of the basis matrix as a chain product of multiple elementary transformation vectors arranged in order. The basis inverse product operation selects the corresponding solution process according to the type of the input vector. When the input is a column vector, a forward solution process is adopted. The scattering rank-one update operation is performed on the input column vector in the order of generation of the elementary transformation vectors to obtain the product of the inverse matrix of the basis matrix and the input column vector. Each step of the scattering rank-one update operation includes: extracting the component value of the input column vector at the pivot position as a scaling factor, clearing the pivot position component to zero, and superimposing the product of the scaling factor and the corresponding elementary transformation vector onto all components of the input column vector. When the input is a row vector, a reverse solution process is adopted. The inner product assignment operation is performed on the input row vector in reverse order of the generation order of the elementary transformation vectors to obtain the product of the input row vector and the inverse of the basis matrix. Each step of the inner product assignment operation includes calculating the inner product value of the input row vector and the current elementary transformation vector and assigning the inner product value to the component of the input row vector at the corresponding pivot position, while the other components of the input row vector remain unchanged.

5. The method according to claim 4, characterized in that, The process of performing parsimony cost calculations on all non-basic variables includes: Through the reverse solution process, corresponding components are extracted from the cost vector of the objective function based on the basic variable index set to form the cost coefficient vector of the current basic variable. The cost coefficient vector is used as the input row vector to obtain the dual variable vector. For each non-basic variable, a parallel execution thread is assigned. Each parallel execution thread reads the non-zero elements of the corresponding constraint matrix column from the high-bandwidth video memory and performs a sparse inner product operation with the dual variable vector. The result of the sparse inner product operation is subtracted from the corresponding original cost coefficient to obtain the simplified cost value of the non-basic variable.

6. The method according to claim 1, characterized in that, The process of obtaining the global minimum parsimony cost value and the corresponding entering base column index adopts a two-level parallel reduction structure, including: In the first-level reduction, the simplified cost values ​​of all non-basic variables obtained by the simplified cost operation are evenly distributed to multiple thread blocks. Multiple threads in each thread block perform tree-structured comparison reduction in shared memory to obtain local minimum values ​​and their local indices. In the second-level reduction, a global reduction is performed on the local minimum values ​​output by all the thread blocks to obtain the global minimum parsimony cost value and the advancing base column index.

7. The method according to claim 1, characterized in that, The extraction of the column vector corresponding to the base column index includes: When the entering basic column index points to the original decision variable column, a parallel execution thread is allocated to each row of the constraint matrix. Each parallel execution thread performs a binary search in the column index range of the compressed sparse row format of the corresponding row to locate the non-zero element value corresponding to the entering basic column index, and assembles the search results of all rows into a dense column vector. When the entering basis column index exceeds the range of the original decision variable columns of the constraint matrix, the column vector corresponding to the entering basis column index is the standard basis vector, and the standard basis vector is directly generated as the dense column vector.

8. The method according to claim 1, characterized in that, The minimum ratio test uses the Brand inverse loop rule. When there are multiple basic variable rows whose ratios are all equal to the minimum ratio, the row with the smallest index value in the basic variable index set is selected as the off-base row index.

9. The method according to claim 4, characterized in that, Also includes: Monitor the cumulative number of elementary transformation vectors in the product form basis inverse representation. When the cumulative number exceeds a preset reconstruction threshold, perform partial principal triangulation on the current basis matrix to regenerate the product form basis inverse representation and reset the cumulative number.

10. The method according to claim 9, characterized in that, Also includes: After each iteration of a preset verification interval, a basic variable column index is selected from the basic variable index set. The basic inverse product operation is performed on the constraint matrix column vector corresponding to the basic variable column index through the forward solution process to obtain the verification vector. The infinite norm residual between the verification vector and the standard basis vector corresponding to the position of the basic variable column index in the basis matrix is ​​compared with a preset precision threshold. When the infinite norm residual exceeds the preset precision threshold, the partial principal component triangular decomposition is triggered.

Images