quicwind/saddlepoint/saddle_sparse.m

   1 % matrix implementation of the saddlepoint method
   2
   3 % repeated subexpressions
   4 Bt=B.';
   5 Dt=D.';
   6 Ct=C.';
   7 ib=inv(Bt*A*B);
   8 ibd=inv(D*ib*Dt);
   9
  10 % solution matrices
  11 P=ib - ib*Dt*ibd*D*ib;
  12 M =C*P*Ct;
  13 rhs=C*P*Bt*A*v0;
  14
  15 % solution components
  16 q=M\rhs;
  17 v=P*(Bt*A*v0-Ct*q);
  18 p=ibd*D*ib*(Bt*A*v0-Ct*q);
  19
  20 % test
  21 res=[Bt*A*B*v + Dt*p + Ct*q - Bt*A*v0
  22     D*v
  23     C*v];
  24 err=norm(res,inf)