binv

(tmp_for_tar/misc/binv.m)

Function Synopsis

b = binv(a,invA)

Help text

       b = binv(a,invA)    block inversion

   b is the inverse of a = [ A  B ],  when invA == inv(A) is known.
                           [ C  D ]   

  This func may be quicker than inv(a), if  A is big, but D is small.

 Last modified: April 2001

Cross-Reference Information

test_binv      tmp_for_tar/misc/test_binv.m

Listing of function file tmp_for_tar/misc/binv.m

##       b = binv(a,invA)    block inversion
##
##   b is the inverse of a = [ A  B ],  when invA == inv(A) is known.
##                           [ C  D ]   
##  
##  This func may be quicker than inv(a), if  A is big, but D is small.
##

## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>
## Last modified: April 2001

function b = binv(a,invA)

P = size(a,1) ;
Q = size(invA,1) ;

#  B = a(1:Q,Q+1:P) ;
#  C = a(Q+1:P,1:Q) ;
#  S = inv( a(Q+1:P,Q+1:P) - C*invA*B) ;
#  B = invA*B ;
#  C = C*invA ;
## i1 = 1:Q ;                   # Maybe a little bit slower
## i2 = Q+1:P ;
B = a(1:Q,Q+1:P) ;
C = a(Q+1:P,1:Q) ;
if (rk = rank(a(Q+1:P,Q+1:P) - C*invA*B)) < P-Q,
   printf("binv : rank is %i rather than %i\n",rk,P-Q) ;
##  keyboard
end
S = pinv( a(Q+1:P,Q+1:P) - C*invA*B, 1e-4) ;

B = invA*B ;			# Don't need to do it, since B goes
				# straight to b. But this seems in fact
				# quicker
C = C*invA ;

## More readable, but slower
## b = [ invA + B*S*C, -B*S ; -S*C, S ] ;

## Less readable, but seems quicker
b = zeros(P) ;
b(1:Q,Q+1:P) = -B*S ;
b(1:Q,1:Q) = invA - b(1:Q,Q+1:P)*C ;
b(Q+1:P,1:Q) = -S*C ;
b(Q+1:P,Q+1:P) = S ;

## B = a(1:Q,Q+1:P) ;
## C = a(Q+1:P,1:Q) ;
#  S = inv( a(Q+1:P,Q+1:P) - a(Q+1:P,1:Q)*invA*a(1:Q,Q+1:P) ) ;
#  a(1:Q,Q+1:P) = invA*a(1:Q,Q+1:P) ;
#  a(Q+1:P,1:Q) = a(Q+1:P,1:Q)*invA ;

#  b = [ invA + a(1:Q,Q+1:P)*S*a(Q+1:P,1:Q), -a(1:Q,Q+1:P)*S ; -S*a(Q+1:P,1:Q), S ] ;