¯
    term = (eye(rows(s))-d0*invpd(d0'*w*d0)*d0'*w);
    /* term is idempotent. Hence  it has all zeros and 1's as eigenvectors.
       its eigenvalue decomposition is ve*lam*inv(ve) so its pseudo
       inverse is ve*lam*inv(ve) = itself. */

inv(term*s*term') = inv(term')*inv(s)*inv(term)
       = term'*inv(s)*term */

    varg = term*s*term';
    {invar,doff1} = pseudoin(varg,-99);
     invar2 = term*invpd(s)*term;
    ¯ output file = trash.out reset;  format /RE 6,2;
     ;output off; ®
     ¯