qqr.src

没有说明

        clear v;
        q = zeros(n,n);     /* initialize space for q */
    else;
        v = seqa(1,1,p);    /* use to create mask */
        v = v .<= v';
        r = x' .* v ;        /* R */
        clear v;
        q = zeros(n,n);     /* initialize space for q */
    endif;

    y = 1~zeros(1,n);


#ifDLLCALL
#else

    if rows(_qrsl) /= 455 or _qrsl[1] $== 0;
        _qrsl = zeros(455,1);
        loadexe _qrsl = qrsl.rex;
    endif;

#endif

    i = 1;
    do until i > k;         /* Compute the Q of QR fame */

#ifDLLCALL
#else

        callexe _qrsl(x,n,n,k,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

#ifDLLCALL

        dllcall qrsl(x,n,n,k,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

        q[.,i] = qy;
        y = shiftr(y,1,0);          /* shift the 1 down 1 row */
        i = i + 1;
    endo;
    retp(q,r,pvt);
endp;

/*
**> qqrep
**
**  Purpose:    Computes the orthogonal-triangular (QR)
**              decomposition of a matrix X, such that:
**
**                          X[.,E] = Q1*R.                      (14)
**
**  Format:     { Q1,R,E } = qqrep(X,PVT);
**
**  Input:      X      NxP matrix.
**
**              PVT    Px1 vector, controls the selection of the pivot
**                     columns:
**
**                          if PVT[i] gt 0 then X[i] is an initial column.
**                          if PVT[i] eq 0 then X[i] is a free column.
**                          if PVT[i] lt 0 then X[i] is a final column.
**
**                     The initial columns are placed at the beginning
**                     of the matrix and the final columns are placed
**                     at the end.  Only the free columns will be moved
**                     during the decomposition.
**
**  Output:     Q1     NxK unitary matrix, K = min(N,P).
**
**              R      LxP upper triangular matrix, L = min(N,P).
**
**              E      Px1 permutation vector.
**
**  Remarks:   Given X[.,E], where E is a permutation vector that permutes
**             the columns of X, there is an orthogonal matrix Q such that
**             Q' * X[.,E] is zero below its diagonal, i.e.,
**
**                   Q'* X[.,E] = [ R ]                           (14a)
**                                [ 0 ]
**
**             where R is upper triangular.
**             If we partition
**
**                   Q = [ Q1 Q2 ]                               (15)
**
**             where Q1 has P columns then
**
**                    X[.,E] = Q1 * R                            (16)
**
**             is the qr decomposition of X[.,E].
**
**             qqrep allows you to control the pivoting.  For example,
**             suppose that X is a data set with a column of ones in the
**             first column.  If there are linear dependencies among the
**             columns of X, the column of ones for the constant may get
**             pivoted away.  This column can be forced to be included
**             among the linearly independent columns using pvt.
**
**             If you want only the R matrix, see qrep.  Not computing Q1
**             can produce significant improvements in computing time and
**             memory usage.
**
**  Globals:    _qrdc, _qrsl
**
**  See Also:   qqr, qre, olsqr
*/

proc (3) = qqrep(x,pvt);
    local flag,n,p,qraux,work,dum,info,job,qy,r,v,q,i,y,k,dif;

    /* check for complex input */
    if iscplx(x);
        if hasimag(x);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            x = real(x);
        endif;
    endif;

    if iscplx(pvt);
        if hasimag(pvt);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            pvt = real(pvt);
        endif;
    endif;

    n = rows(x);
    p = cols(x);
    qraux = zeros(p,1);
    work = qraux;

    dum = 0;
    info = 0;
    job = 10000;    /* compute qy only */
    qy = zeros(n,1);
    x = x';

    flag = 1;

#ifDLLCALL
#else

    if rows(_qrdc) /= 647 or _qrdc[1] $== 0;
        _qrdc = zeros(647,1);
        loadexe _qrdc = qrdc.rex;
    endif;
    callexe _qrdc(x,n,n,p,qraux,pvt,work,flag);

#endif

#ifDLLCALL

    dllcall qrdc(x,n,n,p,qraux,pvt,work,flag);

#endif

    k = minc(n|p);
    dif = abs(n-p);

    if n > p;
        r = trimr(x',0,dif);
        v = seqa(1,1,p);    /* use to create mask */
        r = r .*( v .<= v' );       /* R */
        clear v;
        q = zeros(n,p);     /* initialize space for q */
    elseif p > n;
        v = seqa(1,1,p);    /* use to create mask */
        v = v .<= v';
        v = trimr(v,0,dif);
        r = x' .* v ;        /* R */
        clear v;
        q = zeros(n,n);     /* initialize space for q */
    else;
        v = seqa(1,1,p);    /* use to create mask */
        v = v .<= v';
        r = x' .* v ;        /* R */
        clear v;
        q = zeros(n,n);     /* initialize space for q */
    endif;

    y = 1~zeros(1,n);


#ifDLLCALL
#else

    if rows(_qrsl) /= 455 or _qrsl[1] $== 0;
        _qrsl = zeros(455,1);
        loadexe _qrsl = qrsl.rex;
    endif;

#endif

    i = 1;
    do until i > k;         /* Compute the Q of QR fame */

#ifDLLCALL
#else

        callexe _qrsl(x,n,n,k,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

#ifDLLCALL

        dllcall qrsl(x,n,n,k,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

        q[.,i] = qy;
        y = shiftr(y,1,0);          /* shift the 1 down 1 row */
        i = i + 1;
    endo;
    retp(q,r,pvt);
endp;

/*
**> orth
**
**  Purpose:    Computes an orthonormal basis for the column space
**              of a matrix.
**
**  Format:     y = orth(x);
**
**  Input:      x           NxP matrix.
**
**              _orthtol    global scalar, the tolerance for testing if
**                          diagonal elements are approaching zero.  The
**                          default value is 1.0e-14.
**
**  Output:     y           NxL matrix such that y'y == eye(L) and whose
**                          columns span the same space as the columns of x.
**                          L is the rank of x.
**
**  Globals:    _orthtol, _qrdc, _qrsl
*/

proc orth(x);
    local flag,dif,n,p,qraux,work,rd,pvt,job,karg,dum,info,qy,q,i,y;

    /* check for complex input */
    if iscplx(x);
        if hasimag(x);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            x = real(x);
        endif;
    endif;

    n = rows(x);
    p = cols(x);
    qraux = zeros(p,1);
    work = qraux;
    pvt = qraux;
    flag = 1;       /* use pivoting */

    dum = 0;
    info = 0;
    job = 10000;    /* compute qy only */
    qy = zeros(n,1);

    x = x';

#ifDLLCALL
#else

    if rows(_qrdc) /= 647 or _qrdc[1] $== 0;
        _qrdc = zeros(647,1);
        loadexe _qrdc = qrdc.rex;
    endif;
    callexe _qrdc(x,n,n,p,qraux,pvt,work,flag);

#endif

#ifDLLCALL

    dllcall qrdc(x,n,n,p,qraux,pvt,work,flag);

#endif

    dif = abs(n-p);
    if n > p;
        rd = abs(diag(trimr(x',0,dif)));    /* abs of diagonal of R  */
        karg = sumc( rd .> _orthtol*rd[1,1] );      /* number of diagonal
                                                    :: elements of R greater
                                                    :: than tolerance
                                                    */
        q = zeros(n,p);     /* initialize space for q */
    elseif p > n;
        rd = trimr(x,0,dif);
        rd = abs(diag(rd'));        /* abs of diagonal of R */
        karg = sumc( rd .> _orthtol*rd[1,1] );      /* number of diagonal
                                                    :: elements of R greater
                                                    :: than tolerance
                                                    */
        q = zeros(n,n);     /* initialize space for q */
    else;
        rd = abs(diag(x'));         /* abs of diagonal of R */
        karg = sumc( rd .> _orthtol*rd[1,1] );      /* number of diagonal
                                                    :: elements of R greater
                                                    :: than tolerance
                                                    */
        q = zeros(n,n);     /* initialize space for q1 */
    endif;

    if n == 1;
        y = 1;
    else;
        y = 1~zeros(1,n-1);
    endif;

    i = 1;

#ifDLLCALL
#else

    if rows(_qrsl) /= 455 or _qrsl[1] $== 0;
        _qrsl = zeros(455,1);
        loadexe _qrsl = qrsl.rex;
    endif;

#endif

    do until i > karg;

#ifDLLCALL
#else

        callexe _qrsl(x,n,n,karg,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

#ifDLLCALL

        dllcall qrsl(x,n,n,karg,qraux,y,qy,dum,dum,dum,dum,job,info);

#endif

        q[.,i] = qy;
        y = shiftr(y,1,0);          /* shift the 1 down 1 row */
        i = i + 1;
    endo;
    retp(q[.,1:karg]);
endp;