Posted on behalf of Paras:

Hi,

I am hoping that the c gurus can help me understand how OpenMx matrices

are stored, updated and passed around.

Here is the situation:

We are trying design data-structures for the block-sparse BLAS routines

and are struggling with the issue of ownership of "double" matrices. Our

key data-structure is called the "blockSparseMatrix". For example, there

may be a huge "factor-loading" matrix where each dense sub-block would

map on to one or more OpenMx parameter matrices. In the ordinary

factor-model, each dense sub-block would point to the single

factor-loading matrix. Alternatively in the growth-curve case with

definition variables, each dense-block would point to person specific

factor-loading matrix.

The struct blockSparseMatrix has a field **dns or a "pointer to an array

of doubles representing the dense matrix". From a memory and efficiency

perspective, it makes sense for OpenMx to create, update and destroy the

memory for dense-matrices. Pointer to the dense matrix would be passed

to the backend at the beginning and would be retained throughout the

estimation process as a member of the blockSparseMatrix structure.

>From an OO perspective, having the main application write directly to a

struct-member seems to be problematic. However, I see no alternative to

doing so -- without copying huge amounts of data at each iteration.

>From OpenMx's design perspective and how you envision development of

other plug-ins, what approach would you recommend?

How does OpenMx deal with R? If I understand correctly, R insists that

all data be "copied" between R and external-packages, unless the access

is read-only. Is this what happens at each iteration?

The "Writing R Extensions" manual is a good place to start on learning on memory management interacts with C calls: http://cran.r-project.org/doc/manuals/R-exts.html#Handling-R-objects-in-C. C functions receive arguments as SEXPs, which are pointers to structures with typedef SEXPREC. This leads me to believe that data is not copied when passed over the fence from R to C. Also C programmers needs to explicitly protect data from garbage collection, which suggests that memory is not copied (because it can be garbage collected).

With regards to the decision to keep the front-end as dense matrices, here is my perception of the advantages and disadvantages of this decision. Advantages: it is simpler with respect to interfacing to the front-end. No changes need to be made in R. Disadvantages: there are costs associated with storing dense matrices. There is a time cost associated with converting from dense to sparse at the beginning of optimisation, and sparse to dense at the end of the operation. There is also a memory cost associated with carrying around a huge matrix in the front-end. I believe the memory cost is greater than the time cost. OTOH, I think designing the system in this fashion is good for a first attempt. We should probably be using something akin to the Matrix package in the front-end anyhow, but we're not.

Some design hints on writing a matrix package for the R front-end: (1) make sure your matrices can hold numeric, logical, and character types. The Matrix package can hold numerics and logicals. This handles almost all our cases, except for the labels matrix which holds character vectors, ie. strings. (2) make sure you can specify what is the sparse value. The values matrix holds numeric types where the sparse value is 0. The lbound and ubound matrices hold numeric types where the sparse value is NA.

Will do.

additional questions/thoughts in my reply to Tim's response.

Sparse vs Dense matrices in the front-end:

I think keeping the front-end matrices as dense along with the flexibility/functionality that openMx matrices makese sense. Probably much faster than sparse matrix operations for most common SEM problems.

R insists that its own data remain unmodified by external packages. That holds for (essentially) any memory allocated by R and passed to the external package. Within the external package run (for example, within a call to the C back-end), memory can be allocated, modified, and freed almost without restriction.

If this is a situation where it is preferable to have directly manipulable data, it can be duplicated once at the beginning of the optimization routine, manipulated during each iteration of the optimizer, and then the duplicate destroyed at the end of the routine, all within the back-end. That way the time per iteration isn't increased.

The omxMatrix object at the back-end works this way, using R's memory directly whenever only read access is needed. It is worth noting that the R accessor functions do not guarantee that they will not duplicate data in order to pass it from R to C, though in general they do not.

The whole thing is complicated, of course, by the fact that BLAS tends to operate in-place, and therefore often modifies arguments during computation. As a result, there are many times when duplication of data is necessary at every iteration. When the computations overwrite the original values, though, I'm not sure that there's any way to avoid at least some copying.

That was helpful.

For any direct interface with R, copying will generally be needed, unless the access is read-only. I suppose similar assumptions can be made for the omxModel object passed to the objectiveFunction.

Our sparse-matrix structures and routines will work with pointers to externally allocated dense matrices.

For our backend, we plan to construct our sparse-matrix objects at the first iteration - essentially grabbing pointers to the numeric parts of all omxMatrices. Then at every iteration, we will just need to update pointers to these matrices - if the locations have changed. If no reallocation has occured within OpenMX then our back-end structures will be automatically updated once openMx updates the model object for the objectiveFunction. Either of these options would be better than copying numeric values from omxMatrices to our back-end structures, which is not terribly worse.

Which of these three scenarios should we assume for our back-end?

I would recommend checking/updating your pointers at each iteration.

omxMatrix objects that are passed from the front-end are statically allocated, and will not change their pointer structure, so pointers to them can be used without worrying about changes. omxAlgebras, at present, are not guaranteed to maintain their data pointers, although they should only change if the data matrices grow between iterations, which is not a case I've seen yet.

If you maintain the omxMatrix* pointer to either an omxAlgebra or an omxMatrix, though, you can always access the current location of the data object (or of a specific element in the data object) from that. The function to do so is currently called double* omxLocationOfMatrixElement(omxMatrixObject om, int row, int col), but that is subject to change until the back-end API is finalized.

It's also worth noting that omxObjective objects in the back-end are set up through the call to an initialization function, which is called after all matrices and (hopefully) all dependent algebras have been generated. This is primarily intended to allow objective objects to unpack their arguments, link to their associated data objects, etc., but also provides the opportunity to grab the pointers you need and perform precalculations.

That is what I was hoping. Thanks,.