Applications of QR, RUPDT, and RDODT Subroutines

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: QR                                                  */
 /*   TITLE: Applications of QR, RUPDT, and RDODT Subroutines    */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX                                              */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: MZC                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

proc iml;

 a={  -74   80   18  -11   -4   -8,
       14  -69   21   28    0    7,
       66  -72   -5    7    1    4,
      -12   66  -30  -23    3   -3,
        3    8   -7   -4    1    0,
        4  -12    4    4    0    1};

          /* Second column of b is orthogonal to the columns of a */
          /* third column of b is equal the sum of the first two  */
 b= { 51    -56     -5,
     -61     52     -9,
     -56    764    708,
      69   4096   4165,
      10 -13276 -13266,
     -12   8421   8409};

 aa = a[,1:3];
 m= nrow(aa);
 n= ncol(aa);
 p= ncol(b);

 call QR(q,r,piv,lindep,aa);

 r0 = j(n,n,0.);
 qtb = j(n,p,0.);
 ssq = j(1,p,0.);
 bb=b;
 CALL RUPDT(rup,bup,sup,r0,aa,qtb,bb,ssq);

   /* r and rup should be identical(maybe in different sign) */
 print r rup;

   /* check should be close to bup (maybe in different sign */
 check=(q`*b)[1:3,];
 print check bup;

   /* get ssq from from the last 3 rows of q`*b */
 check=(q`*b)[4:6,];
 check=sqrt(diag(check`*check));

   /* diagonal elements of check should be close to sup */
 print check sup;


   /* downdateing the last row of aa */
 z=aa[6,]; bb=b[6,];
 CALL Rdodt(def,rup1,bup1,sup1,rup,z,bup,bb,sup);

   /* def should be 0 */
 print def;

 aa = a[1:5,1:3];
 m= nrow(aa);
 n= ncol(aa);
 p= ncol(b);

   /* updateing the first 5 rows of aa recursively */
 r = j(n,n,0.);
 qtb = j(n,p,0.);
 ssq = j(1,p,0.);
 do i = 1 to m;
   z = aa[i,1:n]; bb = b[i,1:p];
   CALL RUPDT(rup,bup,sup,r,z,qtb,bb,ssq);
   r = rup; qtb = bup; ssq = sup;
 end;

   /* updateing the first 5 rows of aa at one time */
 r = j(n,n,0.);
 qtb = j(n,p,0.);
 ssq = j(1,p,0.);
 z = aa[1:5,1:n]; bb = b[1:5,1:p];
 CALL RUPDT(rup2,bup2,sup2,r,z,qtb,bb,ssq);


  /* each group of matrices should be identical */
 print sup sup1 sup2, rup rup1 rup2, bup bup1 bup2;

proc iml;

 a={  -74   80   18  -11   -4   -8,
       14  -69   21   28    0    7,
       66  -72   -5    7    1    4,
      -12   66  -30  -23    3   -3,
        3    8   -7   -4    1    0,
        4  -12    4    4    0    1};

        /* Second column of b is orthogonal to the columns of a */
        /* third column of b is equal the sum of the first two  */
 b= { 51    -56     -5,
     -61     52     -9,
     -56    764    708,
      69   4096   4165,
      10 -13276 -13266,
     -12   8421   8409};

 aa = a[1:2,1:3]//a[5:6,1:3];
 m= nrow(aa);
 n= ncol(aa);
 p= ncol(b);

 call QR(q,rr,piv,lindep,aa);
 check=q`*(b[1:2,]//b[5:6,]);
 checkssq=sqrt(check[4,]##2);

 r0 = j(n,n,0.);
 qtb = j(n,p,0.);
 ssq = j(1,p,0.);
 bb=b[1:2,]//b[5:6,];
 CALL RUPDT(rup,bup,sup,r0,aa,qtb,bb,ssq);

   /* downdateing the middle 2 rows of a[,1:3] */
 CALL RUPDT(rupp,bup,supp,r0,a[,1:3],qtb,b,ssq);
 z=a[3:4,1:3]; bb=b[3:4,];
 CALL Rdodt(def,rup1,bup1,sup1,rupp,z,bup,bb,supp);

   /* rr, rup and rup1 should be identical(maybe in different   */
   /*           sign),                                          */
   /* checkssq, sup and sup1 should be same, def should be 0    */
 print rr rup rup1, checkssq sup sup1, def;