Fortran - Using SLATEC BLAS

Author

Using SLATEC BLAS

Daniel Lohmann

2007-06-22, 10:07 pm

I would like to evaluate the quadratic

x'Qx

using SLATEC (or comparable) where x is a vector and Q is a covariance
matrix that is symmetrical and probably positive definite.

All the procedures I've seen for vector-matrix multiplication are similar to
the form: alpha*A*x + .

In general how can I handle a form such as alpha*x*A. assuming I would
proceed as (x'Q)x? It seems that I might transform the vector x into a
matrix X then use something like DSYMM.f but that approach seems wasteful.

andy2O

2007-06-24, 8:09 am

On Jun 22, 3:50 pm, "Daniel Lohmann" <dan...@losoft.de> wrote:
> I would like to evaluate the quadratic
>
> x'Qx
>
> using SLATEC (or comparable) where x is a vector and Q is a covariance
> matrix that is symmetrical and probably positive definite.
>
> All the procedures I've seen for vector-matrix multiplication are similar to
> the form: alpha*A*x + .
>
> In general how can I handle a form such as alpha*x*A. assuming I would
> proceed as (x'Q)x? It seems that I might transform the vector x into a
> matrix X then use something like DSYMM.f but that approach seems wasteful.

the obvious question is why do you want to do (x'Q)x when surely
x'(Qx) is mathematically equivalent? that's just a simple (symmetric)
matrix vector multiply followed by a dot product isn't it, which is
easily accomplished with BLAS?

(x'Q)x = ((x'Q)x)' = x'(x'Q)' = x'(Q'x) = x'(Qx) as Q symmetric

note BLAS does support alpha*A'y for matrix A, vector y and scalar
alpha. mathematically, the result of (alpha*y'*A) is the transpose of
(alpha*A'*y) but as the results are usually stored in 1-D arrays in
fortran, there is no practical difference (you *cannot* tell the
difference between a 1-d row vector and a 1-d column vector in fortran
- they're just 1-d arrays). so blas can do what you want - however as
your Q is symmetric this is unnecessary in your case anyway...

hope that helps,
andy

ps: please note that if performance is important to your code and the
matrices are large you should investigate the ATLAS blas and Goto's
blas (google search will find both) or intel's mkl library which will
provide highly optimised versions of the blas library which may be
much faster than compiling the source code yourself.

Daniel Lohmann

2007-06-24, 8:09 am

"andy2O" <andy2O@hotmail.com> wrote in message
news:1182680647.000057.201660@q69g2000hsb.googlegroups.com...
| On Jun 22, 3:50 pm, "Daniel Lohmann" <dan...@losoft.de> wrote:
| > I would like to evaluate the quadratic
| >
| > x'Qx
| >
| > using SLATEC (or comparable) where x is a vector and Q is a covariance
| > matrix that is symmetrical and probably positive definite.
| >
| > All the procedures I've seen for vector-matrix multiplication are
similar to
| > the form: alpha*A*x + .
| >
| > In general how can I handle a form such as alpha*x*A. assuming I would
| > proceed as (x'Q)x? It seems that I might transform the vector x into a
| > matrix X then use something like DSYMM.f but that approach seems
wasteful.
|
| the obvious question is why do you want to do (x'Q)x when surely
| x'(Qx) is mathematically equivalent? that's just a simple (symmetric)
| matrix vector multiply followed by a dot product isn't it, which is
| easily accomplished with BLAS?
|
| (x'Q)x = ((x'Q)x)' = x'(x'Q)' = x'(Q'x) = x'(Qx) as Q symmetric
|
| note BLAS does support alpha*A'y for matrix A, vector y and scalar
| alpha. mathematically, the result of (alpha*y'*A) is the transpose of
| (alpha*A'*y) but as the results are usually stored in 1-D arrays in
| fortran, there is no practical difference (you *cannot* tell the
| difference between a 1-d row vector and a 1-d column vector in fortran
| - they're just 1-d arrays). so blas can do what you want - however as
| your Q is symmetric this is unnecessary in your case anyway...
|
| hope that helps,
| andy
|
| ps: please note that if performance is important to your code and the
| matrices are large you should investigate the ATLAS blas and Goto's
| blas (google search will find both) or intel's mkl library which will
| provide highly optimised versions of the blas library which may be
| much faster than compiling the source code yourself.
|

Thanks Andy. Your response sure does help a lot.

I am doing x'(Qx) as you pointed by calling DSYMV then DDOT.

I will definitely look into ATLAS, etc. While my cases involve an order of
about 10,000, I will need to do several tens of thousands of these
multiplications.

Sponsored Links