[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: matrix inversion
From: |
Gorazd Brumen |
Subject: |
Re: matrix inversion |
Date: |
Thu, 23 Nov 2006 10:41:27 +0100 |
User-agent: |
Thunderbird 1.5.0.8 (X11/20061116) |
Hi David,
Yes, it is actually qute easy to prove that
the inverse is a series of t NxN full matrix inversions.
Since running dmperm gives me the following result
dmperm (A)
error: dmperm: not available in this version of Octave
(could it be that I am running octave 2.9.7 with forge for (up to)
2.9.6???)
I checked it in matlab though and it is as I have written:
t NxN inversions, you just got confused because in your
example t and N are the same.
I was astonished to there is already a procedure to do what I was
doing by hand. I will use it, it simplifies
my code significantly.
If I get dmperm running in Octave, I would like to try your code, and
once again, thanks for the kind replies.
Regards,
Gorazd
David Bateman wrote:
> Gorazd Brumen wrote:
>> Hello,
>>
>> I have a square matrix that is composed of N^2 block matrices, where
>> every block matrix is relatively large t x t diagonal matrix (both
>> N and t large), i.e.
>>
>> A = [ A_{11}, A_{12} , ... A_{1N}; A_{21} ... A_{2N} ...],
>> where every A_{ij} is a diagonal matrix of size txt.
>>
>> My question is the following: How many operations would
>> spinv need to do this inversion?
>>
>>
>> Thanks and regards,
>> Gorazd
>>
>>
> Ok, just for the hell of it I ran
>
> N=5;
> t=5;
> dt = speye(t);
> A = repmat(dt,N,N);
> [p,q,r,s] = dmperm(A);
> spy(A(p,q))
>
> and I saw that the block triangular factorization of this beast is in
> fact N t-by-t full matrices. So yes a BTF technique really does make
> sense here.. I believe, that something like
>
> [p,q,r,s] = dmperm(A);
> b = speye(rows(A))(p,:);
> A = A (p,q);
> x = zeros(size(A));
> for k = r-1:-1:1
> # Get k-th block
> k1 = r(k);
> k2 = r(k+1) - 1;
>
> # Solve the system
> x(k1:k2,:) = A(k1:k2,K1:k2) \ b(k1:k2,:);
>
> # off-diagonal block back subsitution
> b (1:k1-1,:) = b (1:k1-1,:) - A (1:k1-1, k1:k2) * x (k1:k2,:) ;
>
> endfor
> x(q,:) = x;
>
> might be an efficient manner to solve this type of problem (at least in
> theory).
>
> D.
>
--
Gorazd Brumen
Mail: address@hidden
WWW: http://www.gorazdbrumen.net/
PGP: Key at http://pgp.mit.edu, ID BCC93240