where Q is an orthogonal matrix (i.e transpose(Q)=inverse(Q)) and R is an upper traingular matrix. then the product QR is called a
* 다음은 Octave 실행 창에서 help 명령으로 qr 사용법을 알아본 것이다.
/***********************************************
Filename: example-05.c
Purpose:
Test GSL library.
http://www.gnu.org/software/gsl/
On Window XP & cygwin:
On Ubuntu
Compile: gcc -Wall -I/usr/local/include -c example-05.c
Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm
On Mac OS X Lion:
Compile: gcc -Wall -I/usr/local/include -c example-05.c
Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm
Execute: ./example-05
Result:
A is a 2 by 2 matrix
A =
[ [ 3.00000000 4.00000000 ],
[ 1.00000000 2.00000000 ] ]
A =
[ [ -3.16227766 -4.42718872 ],
[ 0.16227766 0.63245553 ] ]
tau = [ 1.94868330 0.00000000 ]
R =
[ [ -3.16227766 -4.42718872 ],
[ 0.00000000 0.63245553 ] ]
Q =
[ [ 0.00000000 0.00000000 ],
[ 0.16227766 0.00000000 ] ]
***************************************************/
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
void printVector(char *msg, gsl_vector *x, const char *fmt) {
if (msg != NULL && strlen(msg) > 0)
printf("%s", msg);
size_t n = x->size;
size_t i;
printf("[ ");
printf (fmt, gsl_vector_get(x, 0));
for (i = 1; i < n; i++) {
printf (fmt, gsl_vector_get(x, i));
}
printf(" ]\n");
}
void printMatrix(char *msg, gsl_matrix *a, char *fmt) {
if (msg != NULL && strlen(msg) > 0)
printf("%s\n", msg);
size_t m = a->size1, n = a->size2;
size_t i, j;
printf("[ [ ");
printf (fmt, gsl_matrix_get(a, 0, 0));
for (j = 1; j < n; j++) {
printf (fmt, gsl_matrix_get(a, 0, j));
}
printf(" ]");
if (m > 1)
printf(",\n");
for (i = 1; i < m; i++) {
printf(" [ ");
printf (fmt, gsl_matrix_get(a, i, 0));
for (j = 1; j < n; j++) {
printf (fmt, gsl_matrix_get(a, i, j));
}
printf(" ]");
if (i < m -1)
printf(",\n");
}
printf(" ]\n");
}
void product(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
size_t m = a->size1, n = a->size2, p = b->size2;
size_t i, j, k;
double sum;
for (i = 0; i < m; i++) {
for (k = 0; k < p; k++) {
sum = 0.0;
for (j = 0; j < n; j++) {
sum = sum + gsl_matrix_get(a, i, j)*gsl_matrix_get(b, j, k);
}
gsl_matrix_set (c, i, k, sum);
}
}
}
void productWIthRightVector(gsl_vector *c, const gsl_matrix *a, const gsl_vector *b) {
size_t m = a->size1, n = b->size;
size_t i, j;
double sum;
for (i = 0; i < m; i++) {
sum = 0.0;
for (j = 0; j < n; j++) {
sum = sum + gsl_matrix_get(a, i, j)*gsl_vector_get(b, j);
}
gsl_vector_set(c, i, sum);
}
}
void productWIthLeftVector(gsl_vector *c, const gsl_vector *a, const gsl_matrix *b) {
size_t m = a->size, n = b->size2;
size_t i, j;
double sum;
for (j = 0; j < n; j++) {
sum = 0.0;
for (i = 0; i < m; i++) {
sum = sum + gsl_vector_get(a, j)*gsl_matrix_get(b, i, j);
}
gsl_vector_set(c, j, sum);
}
}
void add(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
size_t m = a->size1, n = a->size2;
size_t i, j;
double v;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
v = gsl_matrix_get (a, i, j) + gsl_matrix_get (b, i, j);
gsl_matrix_set (c, i, j, v);
}
}
}
void subtract(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
size_t m = a->size1, n = a->size2;
size_t i, j;
double v;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
v = gsl_matrix_get (a, i, j) - gsl_matrix_get (b, i, j);
gsl_matrix_set (c, i, j, v);
}
}
}
void multiply(gsl_matrix *c, double x, const gsl_matrix *a) {
size_t m = a->size1, n = a->size2;
size_t i, j;
double v;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
v = x*gsl_matrix_get (a, i, j);
gsl_matrix_set (c, i, j, v);
}
}
}
void set(gsl_matrix *c, const gsl_matrix *a) {
gsl_matrix_memcpy(c, a);
}
double trace(const gsl_matrix *a) {
size_t m = a->size1;
size_t i;
double sum = 0;
for (i = 0; i < m; i++) {
sum = sum + gsl_matrix_get (a, i, i);
}
return sum;
}
void transpose(gsl_matrix *c, const gsl_matrix *a) {
size_t m = a->size1, n = a->size2;
size_t i, j;
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++) {
gsl_matrix_set (c, i, j, gsl_matrix_get (a, j, i));
}
}
}
void getQR(gsl_matrix *Q, gsl_matrix *R, const gsl_matrix *A, const gsl_vector *tau) {
size_t m = A->size1;
size_t n = A->size2;
size_t i, j;
for (i =0 ; i < m; i++) {
for (j =0 ; j < n; j++) {
if (i > j) {
gsl_matrix_set(R, i, j, 0);
gsl_matrix_set(Q, i, j, gsl_matrix_get(A, i, j));
}
else {
gsl_matrix_set(R, i, j, gsl_matrix_get(A, i, j));
gsl_matrix_set(Q, i, j, 0);
}
}
}
}
int main(void) {
size_t n = 2;
gsl_vector *tau = gsl_vector_alloc(n);
// size_t i, j;
gsl_matrix *A = gsl_matrix_alloc(n, n);
gsl_matrix *Q = gsl_matrix_alloc(n, n);
gsl_matrix *R = gsl_matrix_alloc(n, n);
gsl_matrix_set (A, 0, 0, 3);
gsl_matrix_set (A, 0, 1, 4);
gsl_matrix_set (A, 1, 0, 1);
gsl_matrix_set (A, 1, 1, 2);
printf ("A is a 2 by 2 matrix\n");
printMatrix("A = ", A, "%12.8f");
gsl_linalg_QR_decomp(A, tau);
printMatrix("A = ", A, "%12.8f");
printVector("tau = ", tau, "%12.8f");
getQR(Q, R, A, tau);
printMatrix("R = ", R, "%12.8f");
printMatrix("Q = ", Q, "%12.8f");
gsl_matrix_free (Q);
gsl_matrix_free (R);
gsl_matrix_free (A);
gsl_vector_free (tau);
return 0;
}