Re: Transpose(A)*P*A [message #62837 is a reply to message #62831] |
Fri, 10 October 2008 08:03  |
Vince Hradil
Messages: 574 Registered: December 1999
|
Senior Member |
|
|
On Oct 10, 9:25 am, "mapper...@gmail.com" <zjwan...@gmail.com> wrote:
> hello,
>
> I have a question about how to improve the computation speed when deal
> with non-linear equation
> A x = l, P is the weight for each
> row, P is M*M
> M*N N*1 M*1
>
> then I have to build normal matrix which is
> Transpose(A)*P*A x = Transpose(A)*P*l
> N*N N*1 N*1
>
> then x can be solved.
>
> When M is bigger as 20000, N as 3000, the time to build AtPA is almost
> one hour, that is too long. the code is:
>
> PRO NormalMatrix, A, l, ATPA, ATPL, Weight = P
> szA = SIZE(A, /DIMENSIONS)
> if N_ELEMENTS(szA) ne 2 then begin
> mess = dialog_message('Input matrix must be M*N format!', /
> information)
> return
> endif
> szL = SIZE(L, /DIMENSIONS)
> if (N_ELEMENTS(szL) ne 2 and szL[0] ne 1) then begin
> mess = dialog_message('Input l must be 1*M format!', /information)
> return
> endif
> if (szA[1] ne szL[1]) then begin
> mess = dialog_message('Input A and l must be same rows!', /
> information)
> return
> endif
> if not keyword_set(P) then P = fltarr(szA[1])*0.0+1.0
> ATPA = fltarr([szA[0], szA[0]])
> ATPL = fltarr(1, szA[0])
>
> for i=0L, szA[0]-1 do begin
> t1 = systime(1)
> ATPA[i, i:szA[0]-1] = (p*A[i, *])##A[i:szA[0]-1, *]
> ATPA[i:szA[0]-1, i] = ATPA[i, i:szA[0]-1]
> ATPL[0, i] = (p*A[i, *])##l
> t2 = systime(1)
> print, sza[0], i, t2-t1
> endfor
>
> END
>
> I like to know how I can improve it to few seconds.
I think I'm confused - won't be the first time - why not just do the
multiplication?
atpa = matrix_multiply(a,p##a,/atranspose)
or maybe it's
atpa = matrix_multiply(p##a,a,/btranspose)
|
|
|