Linear Algebra and the C Language/a0lk

__NOTOC__ Install and compile this file in your working directory. /* ------------------------------------ */ /* Save as : c03d.c */ /* ------------------------------------ */ #include "v_a.h" /* ------------------------------------ */ /* ------------------------------------ */ #define RB R3 #define CB C1 /* ------------------------------------ */ /* ------------------------------------ */ int main(void) { double b[RB*(CB)]={ +0.333333333333, +0.666666666667, +1.000000000000 }; double x[RB*C1]={ -1, 2, -3 }; double **B = ca_A_mR(b,i_mR(RB,CB)); double **BT = i_mR(CB,RB); double **BTB = i_mR(CB,CB); double **invBTB = i_mR(CB,CB); double **invBTB_BT = i_mR(CB,RB); double **B_invBTB_BT = i_mR(RB,RB); double **Id = eye_mR(i_mR(RB,RB)); double **V = i_mR(RB,RB); /* Id-(B inv(BT B) BT) */ double **X = ca_A_mR(x,i_mR(RB,C1)); double **VX = i_mR(RB,C1); clrscrn(); printf(" B is a basis for the orthogonal complement of AT : \n\n" " Find a transformation matrix for \n" " a projection onto R%d : \n\n" " Proj(x) = [Id-(B inv(BTB)BT)] \n\n" " B: See the previous c file",RB); p_mR(B,S5,P4,C7); stop(); clrscrn(); printf(" BT:"); p_mR(transpose_mR(B,BT),S5,P4,C7); printf(" BT B:"); p_mR(mul_mR(BT,B,BTB),S5,P4,C7); printf(" inv(BTB):"); invBTB[R1][C1] = 1./BTB[R1][C1]; p_mR(invBTB,S5,P4,C7); printf(" inv(BTB) BT :"); p_mR(mul_mR(invBTB,BT,invBTB_BT),S5,P4,C7); printf(" B inv(BTB)BT:"); p_mR(mul_mR(B,invBTB_BT,B_invBTB_BT),S5,P4,C7); printf(" V = Id-(B inv(BTB)BT)"); p_mR(sub_mR(Id,B_invBTB_BT,V),S5,P4,C7); stop(); clrscrn(); printf(" V is transformation matrix for \n" " a projection onto a R%d subspace \n\n" " V:",RB); p_mR(V,S5,P4,C7); printf(" X:"); p_mR(X,S5,P4,C7); printf(" Proj(x) = [Id-(B inv(BTB)BT)] x \n\n" " Proj(x) = V x Verify the result with the first C file"); p_mR(mul_mR(V,X,VX),S5,P4,C7); stop(); f_mR(B); f_mR(BT); f_mR(BTB); f_mR(invBTB); f_mR(invBTB_BT); f_mR(V); /* Id-(B inv(BT B) BT) */ f_mR(X); f_mR(VX); return 0; } /* ------------------------------------ */ /* ------------------------------------ */ Screen output example: B is a basis for the orthogonal complement of AT : Find a transformation matrix for a projection onto R3 : Proj(x) = [Id-(B inv(BTB)BT)] B: See the previous c file +0.3333 +0.6667 +1.0000 Press return to continue. BT: +0.3333 +0.6667 +1.0000 BT B: +1.5556 inv(BTB): +0.6429 inv(BTB) BT : +0.2143 +0.4286 +0.6429 B inv(BTB)BT: +0.0714 +0.1429 +0.2143 +0.1429 +0.2857 +0.4286 +0.2143 +0.4286 +0.6429 V = Id-(B inv(BTB)BT) +0.9286 -0.1429 -0.2143 -0.1429 +0.7143 -0.4286 -0.2143 -0.4286 +0.3571 Press return to continue. V is transformation matrix for a projection onto a R3 subspace V: +0.9286 -0.1429 -0.2143 -0.1429 +0.7143 -0.4286 -0.2143 -0.4286 +0.3571 X: -1.0000 +2.0000 -3.0000 Proj(x) = [Id-(B inv(BTB)BT)] x Proj(x) = V x Verify the result with the first C file -0.5714 +2.8571 -1.7143 Press return to continue.