Newton Excel Bach, not (just) an Excel Blog

Previous Python Post

The analysis of structural frames requires the solution of large sparse matrix systems, which rapidly becomes very time consuming using a VBA based solver, especially for 3D analysis.  I have previously presented versions of my frame analysis spreadsheet that link to compiled solvers (in Fortran, C++, and C#), but these have still required the generation of the complete matrix, consisting mostly of zeros, which needlessly limits the size of model that can be handled.

The Python Pysparse library contains a number of formats for handling and storing sparse matrices, and importantly for frame analysis, an iterative solver for sparse symmetric matrices.   I am now working on adding this solver (and also a sparse direct solver) to my frame analysis spreadsheet.  This post looks at how to generate the matrix data in Excel, and how to link this data to the Python routines, using Pyxll.

The Pysparse iterative solvers are found in the itsolvers module, which includes the following methods: PCG, MINRES, QMRS, BICGSTAB, and CGS.  Some brief experimentation found little difference in performance, but the QMRS was the fastest in my tests, so that is what I have used.  The procedure for using this routine is:

• Import the data from Excel as a Numpy Array.
• Process  the data in linked list format.
• Convert to SSS format.
• Run pre-conditioner routine.
• Run iterative solver routine
• Return the results to Excel as a Numpy_column array.

The code to do all this is:

@xl_func("numpy_array K, numpy_row b, float tol, int its: numpy_column")
def py_ItSolve2(K, b, tol, its):
    n = len(b)
    A = (ll_mat_from_array(K, n)).to_sss()
    x = zeros(n)
    D = precon.ssor(A, 1.0, 1)
    info, iter, relres = itsolvers.qmrs(A, b, x, tol, its, D)
    return x

def ll_mat_from_array(K,n):
    nr = len(K)
    A = spmatrix.ll_mat_sym(n,n)
    for row in range(0,nr):
        i = int(K[row][0])-1
        j = int(K[row][1])-1
        x = K[row][2]
        A[i,j] = A[i,j]+x
    return A

The input for the py_ItSolve2 routine consists of:

• K, the frame stiffness matrix (see below for more details of the required format)
• b, the load vector, a single column with one row for each frame node freedom
• tol, the target tolerance required for the iterative solution
• its, the maximum number of iterations

The stiffness matrix  is defined in three columns, two defining a pair of freedoms (the row and column number of the full matrix), and the third the associated stiffness value.  This data is generated by a VBA routine from the spreadsheet by calculating the freedom values for each beam in turn.  This results in many freedom pairs being repeated, with stiffness values that must be summed to find the resultant stiffness value for the global matrix.  This is done in the ll_mat_from_array routine, which creates a linked list matrix, then simply steps through the input data and adds the stiffness values for each of the freedom pairs.  The linked list array is then converted to sss (Sparse Skyline) format with the method .to_sss().  Finally this array is preconditioned, solved, and the resulting array of node deflections is returned to the calling VBA routine.

Initially I generated the K array as a VBA double array, but discovered that passing this to the Python routines using application.run had a maximum limit of 65536 (2^16) rows.  I experimented with generating the array in Python, but discovered that (surprisingly) this was substantially slower than the same operation in VBA.  Finally I returned to generating the array in VBA, then writing the data to a named spreadsheet range, and passing this range to the Python routine.  The VBA code for this operation is:

Range("kv_").ClearContents
    Range("kv_").Resize(lastk, 3).Name = "kv_"
    Range("kv_").Value = kv

        ' !------------------------equation solution -----------------------------------
    Select Case SolverType
    Case 1
        DefA = Application.Run("py_ItSolve2", Range("kv_"), Loads, Py_Tol, Py_Its)

    Case 2
        DefA = Application.Run("py_LUSolve2", Range("kv_"), Loads)
    End Select

The top of the kv_ range looks like:

kv_ array

A very similar approach is used to run the Pysparse direct sparse solver:

@xl_func("numpy_array K, numpy_row b: numpy_column")
def py_LUSolve2(K, b):
    n = len(b)
    A = (ll_mat_from_array(K, n)).to_csr()
    x = zeros(n)
    LU = superlu.factorize(A, diag_pivot_thresh=0.0)
    LU.solve(b, x)
    return x

In this case the linked list array is converted to a csr (Compressed Sparse Row) format, and the resulting array is then factorised and solved using the superlu routine.