rank-k updating, with packed storage

[chxyike]

Hello all,

I have a question about rank-k updating. Making use of the subroutine 'dsyrk', one can update the symmetric matrix with only half of the matrix is updated. But with this subroutine, the other half part, which is not referenced, is still in the memory, of the case which I want to avoid.

My question is, if I want to update a matrix like C=alfa*A'*A + C, and C is a m-by-m matrix in packed symmetric storage, A is a n by m matrix, how can I implement this operation with BLAS or Lapack library, with high performance.

I know that one can use 'dspr' to realize this updating, but since the Blas level 1 operation is less efficient than level 3, I hope there are some ways to improve the computation speed.

Thanks a lot.

[Julien Langou]

What about giving a try to DSFRK? The data structure for the symmetric matrix is a little tricky so at first you might want to use DTRTTF to convert from full to rectangular full packed (RFP), use DSFRK, then convert back to full with DTFTTR. Please have a look at Section 9 of http://www.netlib.org/lapack/lapack-3.2.html and maybe: Fred G. Gustavson, Jerzy Wasniewski, Julien langou, and Jack J. Dongarra. Rectangular Full Packed Format for Cholesky’s Algorithm: Factorization, Solution and Inversion. LAPACK Working Note 199, April 2008.

We would be happy if you report any success, failure or problem with the routine. The routine is a prototype so we might change the interface based on various needs and feedback. (Yours!) A priori this is the routine you are looking for. Brand new.

Best wishes,
--julien.

[chxyike]

Thanks a lot for your reply.

I will follow your suggestion, using dsfrk for this operation.

Best regards,

[Julien Langou]

Using RFP format should give you Level 3 BLAS performance with half storage. Have a look, let us know if something does not work, etc.

[chxyike]

However, I cannot get the correct result with DSFRK. Would you please give me some example how to call this subroutine?

The A matrix is not symmetric, How can one save it in memory in RFP storage. Please give me some example. Thanks in advance.

[Julien Langou]

Thanks for considering the subroutine.
In DSFRK, A is a "regular" rectangular matrix with leading dimension LDA.
It has the exact same form/shape/data layout /role has the A in DSYRK.
What changes from DSYRK to DSFRK is the matrix C.
In DSFRK, the matrix C is an array of size N*(N+1)/2.

Code: Select all

% A pseudocode to form the normal equations would be something like:
% initialize the M by N matrix A (with leading dimension LDA) as for DSYRK
% allocate the N*(N+1)/2 array C. (It will hold the representation of A^H * A in RFP format)
      CALL  ZHFRK( 'N', 'U', 'C', N, M, DCMPLX(1.0D00,1.0D00), A, LDA, DCMPLX(0.0D00,0.0D00), C )
% Note the first two arguments 'N' and 'U' indicate how you want the RFP matrix stored, there
% are four different possibilities. In a first use of the subroutine, feel free to use any of the four
% combinations.


Note that once you have the normal equations in C in RFP format there are some subroutines to help you move one from there.
You can do the Cholesky factorization of the normal equations for example with zpftrf.f. (While leaving them in RFP format.)
You can convert back to regular packed format zpfttp.f and general full format with zpfttr.f.

[chxyike]

Thanks a lot. Your suggestions is really helpful.

My problem is, if I use zpfttp to change the matrix from RFP to normal packed format, then I need a extra array to the the matrix in normal packed format. Can we save this matrix? Considering the C matrix can reach 9 Gigabytes or even more.

Best regards,

[Julien Langou]

Your problem makes sense to me. That's a problem.
What do you want to do with this 9Go packed matrix?
The idea would be to try to do it with the RFP matrix directly. That should be feasible.
( I am scare of reading of what you want to do !!! )
--julien

(One alternative option is always to start thinking about inplace convertion from TP to TF ... No idea if this is possible. I think so given the recent work of Fred Gustavson.)

[chxyike]

I am processing GOCE satellite data, which dedicating to estimate up to 63000 parameters ( 44 100 at the moment ), with more than 20 736 000 observations.

It would be easy for me to handle the matrix in general packed format, since I need to do some some operation such as pre-elimination.

[Julien Langou]

Hey pretty cool application, you've got there.

A little effort is needed. Can you be more concrete in the stack of functionality you need, please?

So

(1) given an m-by-n matrix of new observations A update the n-by-n normal equations C,

That I think I understand.
(Is the word "update" appropriate?)

(2) pre-elimate

what does this mean? you select a subset of three integers say [1, 10, n] and then you want to "pre-eliminate" so you want the (n-3)-by-(n-3) SPD matrix that is the Schur complement without these variables. Do you know them in advance? Are they last? first?

(3) ??

At some point you want to Cholesky right?

Well all this is very interesting for us. We were also interested in developing a QR like factorizaton with an updated QR routine but that work on a packed R.
(More stable.) It would be nice if we can understand each other on what you need.

The data structure for RFP is not that bad. In some sense it is easier than PP.
See:
http://www-math.cudenver.edu/~langou/la ... 2/rfp.html

Julien.

[chxyike]

Thanks for your kindly asking.

(1). Since the A matrix is huge, I can only use several thousand rows at one time. So the C matrix is updated as
C = 0
C = C + A1'*A1

C = C + A2' * A2

C = C + A3' *A3
......
where A1, A2,A3 ... are sequential rows of the huge matrix A, for example, A has 10 000 rows, A1 could have 1000 rows.
so we call this is updating.

(2) Pre-Eliminate

I am going generate a normal equation every day, then combine all the normal equation in one month or one year together. There are some parameters which is only related to the specific day, and we do not need these parameters. So we pre-Eliminate these parameters, only left the parameters which have the same physical meaning to all the normal equations. Otherwise the scale of the problem will be even larger, which is not necessary.

(3) Cholesky

After I combine ( merge ) all the daily normal equation together, I will do Cholesky factorization, and even compute the inverse of the normal matrix.

Comparing dgemm and dsfrk, which one will be faster?

Best regards from Munich,