Re: [Lapackpp-devel] Cholesky Decomposition A = D'D

[Christian S. <sti...@tu...>]

Am Donnerstag, 24. M=E4rz 2005 22:18 schrieb Jacob (Jack) Gryn:
> > -----Original Message-----
> > Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User
> > Guide", but that's about it.
>
> Strange, I did a google on
> LaSpdFactDouble and
> LaSymmBandFactDouble
>
> It turned up 0 results.

My point is that the only additional (and old) documentation is the one on=
=20
http://math.nist.gov/lapack++/ . Nothing else exists.

> > > As it stands, I don't really see anything to convert from a symmetric
> > > LaGenMatDouble to a LaSymmBandMatDouble or vice versa.
> >
> > A general conversion doesn't make sense, does it? I guess you mean if
> > and only if a general matrix happens to be a symmetric matrix, then you
> > would like to be able to convert it to the LaSymmBandMatDouble?=20
>
> Well, basically, I multiply a few LaGenMatDouble matrices together, and by
> construction the result is a Symmetric, Positive Definite matrix.
>
> Now, this matrix, I need to do a Cholesky factorization on.

Sure, I understand. Of course Lapack[pp] has no means of knowing that this=
=20
result happens to be symmetric, and it probably isn't even exactly symmetri=
c=20
because of rounding errors (on the order of +-1e-15). So in any case you ha=
ve=20
to construct a new LaSymmMatDouble class which you need to tell which part =
of=20
your general matrix should be used. By the way, maybe the constructor which=
=20
takes a (double*), i.e. a pointer to an array of doubles, is the right choi=
ce=20
for this, and A(0,0) will return the pointer to the full array of the=20
original matrix:

LaGenMatDouble A(10,20);
// do something with A; then=20
LaSymmMatDouble B( A(0,0), A.size(0), A.size(1) );

but I'm not sure which part of A is used in B and which one is ignored (upp=
er=20
vs. lower triangular matrix). You would have to look this up in the=20
documentation of the fortran lapack functions that are used.

Again, why are you also talking about the banded matrices?

> Just to test, I created a function in linslv.cc / laslv.h called
> void LaMatFactorize(const LaGenMatDouble &A,LaGenMatDouble& AF)
> based on the function in spdfd.h, but used it with a LaGenMatDouble as
> input.
>
> The result seems to be 'almost' correct.  The function will only look at
> the upper triangle portion of the input matrix (it assumes symmetry), and
> the result is to be an upper-triangular matrix.
>
> The upper triangular portion of the resulting matrix is correct and does
> compare to the results obtained in matlab.  The bottom triangle needs to =
be
> zeroed out.

That results sounds 'exactly' correct to me, because that's the way how=20
fortran lapack stores symmetric matrices internally. One half of them is=20
simply disregarded, so it doesn't even need to be zeroed out. I guess you=20
refer to what you retrieve when you print it to an std::ostream? In that=20
case, the respective operator would have to treat this similar as the lapac=
k=20
functions, i.e. ignoring one half and printing the full matrix solely based=
=20
on the other half.

> So, as a result, I have a working function, the problem is, if I manually
> zero out the bottom triangle, it may not be as optimal as lapack was
> intended to be.
>
> My questions now would be
> a) Is there a quick way of doing this zero-out?

Again, that part is simply ignored. It shouldn't even be zeroed out, simply=
=20
ignored.

> b) Should I submit this as a patch?

No. See above. If you need this special case for yourself, then that's clea=
ry=20
well placed in your own code.

Thanks for asking, anyway. It's nice to see that the library is actually be=
ing=20
used :-))

Christian