## 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;

```