Linear algebra

-->H=[1.    1.    1.    0.;
-->   2.   -1.    0.    1;
-->   1.    0.    1.    1.;
-->   0.    1.    2.   -1];
-->ww=spec(H)
 ww  =
!   2.7320508 !
! - 2.7320508 !
!   0.7320508 !
! - 0.7320508 !
-->//             STABLE SUBSPACES
-->[X,d]=schur(H,'cont');
-->X'*H*X
 ans  =
! - 2.7320508  - 1.110D-15    0.           1.        !
!   0.         - 0.7320508  - 1.         - 7.772D-16 !
!   7.216D-16    0.           2.7320508    0.        !
!   0.         - 6.106D-16    0.           0.7320508 !
-->[X,d]=schur(H,'disc');
-->X'*H*X
 ans  =
!   0.7320508    0.           7.772D-16    1.        !
!   0.         - 0.7320508  - 1.           8.604D-16 !
!   0.           0.           2.7320508  - 1.166D-15 !
!   7.772D-16    1.110D-15  - 1.277D-15  - 2.7320508 !
-->//Selection of user-defined eigenvalues (# 3 and 4 here);
-->deff('[flg]=sel(x)',...
-->'flg=0,ev=x(2)/x(3),...
--> if abs(ev-ww(3))<0.0001|abs(ev-ww(4))<0.00001 then flg=1,end')
-->[X,d]=schur(H,sel)
 d  =
    2.
 X  =
! - 0.5705632  - 0.2430494  - 0.6640233  - 0.4176813 !
! - 0.4176813    0.6640233  - 0.2430494    0.5705632 !
!   0.5705632  - 0.2430494  - 0.6640233    0.4176813 !
!   0.4176813    0.6640233  - 0.2430494  - 0.5705632 !
-->X'*H*X
 ans  =
!   0.7320508    0.           7.772D-16    1.        !
!   0.         - 0.7320508  - 1.           8.604D-16 !
!   0.           0.           2.7320508  - 1.166D-15 !
!   7.772D-16    1.110D-15  - 1.277D-15  - 2.7320508 !
-->//               With matrix pencil
-->[X,d]=gschur(H,eye(H),sel)
 d  =
    2.
 X  =
!   0.5705632    0.2430494    0.6640233    0.4176813 !
!   0.4176813  - 0.6640233    0.2430494  - 0.5705632 !
! - 0.5705632    0.2430494    0.6640233  - 0.4176813 !
! - 0.4176813  - 0.6640233    0.2430494    0.5705632 !
-->X'*H*X
 ans  =
!   0.7320508    0.           9.576D-16    1.        !
!   0.         - 0.7320508  - 1.           0.        !
!   8.882D-16    0.           2.7320508    0.        !
!   0.           0.           0.         - 2.7320508 !
-->//            block diagonalization
-->[ab,x,bs]=bdiag(H);
-->inv(x)*H*x
 ans  =
!   2.7320508    1.610D-15    0.           0.        !
! - 3.664D-15  - 2.7320508    0.           6.661D-16 !
!   0.           0.           0.7320508  - 7.910D-16 !
!   0.           0.           0.         - 0.7320508 !
-->//                     Matrix pencils
-->E=rand(3,2)*rand(2,3);
-->A=rand(3,2)*rand(2,3);
-->s=poly(0,'s');
-->w=det(s*E-A)   //determinant
 w  =
                           2
  - 0.0006293s - 0.0146503s
-->[al,be]=gspec(A,E);
-->al./(be+%eps*ones(be))
 ans  =
!   1.045D+15 !
! - 0.0429579 !
! - 5.377D-17 !
-->roots(w)
 ans  =
!   0         !
! - 0.0429579 !
-->[Ns,d]=coffg(s*E-A);    //inverse of polynomial matrix;
-->clean(Ns/d*(s*E-A))
 ans  =
!   1     0     0  !
!   -     -     -  !
!   1     1     1  !
!                  !
!   0     1     0  !
!   -     -     -  !
!   1     1     1  !
!                  !
!   0     0     1  !
!   -     -     -  !
!   1     1     1  !
-->[Q,M,i1]=pencan(E,A);   // Canonical form;
    rank A^k    rcond
        2.      0.908E+00
    rank A^k    rcond
        2.      0.209E-01
-->clean(M*E*Q)
 ans  =
!   1.    0.    0. !
!   0.    1.    0. !
!   0.    0.    0. !
-->clean(M*A*Q)
 ans  =
!   0.0307333  - 0.0560490    0. !
!   0.0404070  - 0.0736912    0. !
!   0.           0.           1. !