How to multiply a matrix by its transpose [message #190] |
Thu, 19 September 1991 20:54 |
landsman
Messages: 93 Registered: August 1991
|
Member |
|
|
I am trying to determine the best way in IDL to multiply a matrix
by its transpose. The elegant IDL solution, for a matrix A
C = transpose(A) # A ;Elegant IDL code
requires approximately a factor of two more multiplications than are
really needed (because the output is symmetric). This is the
bottleneck in a linear least squares program I am running which
typically requires an hour of execution time. (The same bottleneck
occurs in the User Library procedure CURVEFIT.)
On the other hand, explicitly taking the inner products of the
the needed columns and rows, i.e.
sz = size(A) ;Ugly "Fortran" like code
nterm = sz(2)
c = fltarr(nterm,nterm,/NOZERO )
for k=0,nterm-1 do begin ;Create normal matrix C
for l=0,k do c(k,l) = total(A(*,k)*A(*,l))
endfor
c = c + transpose(c)
is inefficient in IDL because it requires loops over both rows and
columns in the output matrix. Does anyone know of an efficient way
to do this IDL? (I suspect something could be done with the IDL TRED2
routine from "Numerical Recipes" using the "QR decomposition" but I'm
not real familiar with this.)
Thanks,
Wayne Landsman landsman@stars.gsfc.nasa.gov
ST Systems Co.
Code 681
NASA/GSFC
Greenbelt, MD 20771
(301) - 286 - 3625
|
|
|