test/test_qr.c

/*
 * Library:  lmfit (Levenberg-Marquardt least squares fitting)
 *
 * File:     demo/curve1.c
 *
 * Contents: Example for curve fitting with lmcurve():
 *           fit a data set y(x) by a curve f(x;p).
 *
 * Note:     Any modification of this example should be copied to
 *           the manual page source lmcurve.pod and to the wiki.
 *
 * Author:   Joachim Wuttke <j.wuttke@fz-juelich.de> 2004-2013
 *
 * Licence:  see ../COPYING (FreeBSD)
 *
 * Homepage: apps.jcns.fz-juelich.de/lmfit
 */

#include "lmcurve.h"
#include <stdio.h>
#include <stdlib.h>

void lm_qrfac(int m, int n, double *a, int *ipvt,
              double *rdiag, double *acnorm, double *wa);

void set_orthogonal( int n, double *Q, double* v )
{
    int i, j;
    double temp = 0;
    for (i=0; i<n; ++i)
        temp += v[i]*v[i];
    for (i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            Q[j*n+i] = ( i==j ? 1 : 0 ) - 2*v[i]*v[j]/temp;
}

void matrix_multiplication( int n, double *A, double *Q, double *R )
{
    int i,j,k;
    double temp;
    for (i=0; i<n; ++i) {
        for (j=0; j<n; ++j) {
            temp = 0;
            for (k=0; k<n; ++k) {
                temp += Q[k*n+i]*R[j*n+k];
            }
            A[j*n+i] = temp;
        }
    }
}

void matrix_transposition( int n, double *A )
{
    int i,j;
    double temp;
    for (i=0; i<n; ++i) {
        for (j=i+1; j<n; ++j) {
            temp = A[j*n+i];
            A[j*n+i] = A[i*n+j];
            A[i*n+j] = temp;
        }
    }
}

void print_matrix(int m, int n, double *a)
{
    int i,j;
    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            printf( "%8.4g ", a[j*m+i] );
        }
        printf ("\n");
    }
}

int main( int argc, char *argv[] )
{
    double A[9], rdiag[3], acnorm[3], wa[3];
    int i, ipvt[3];

    if ( argc!= 10 ) {
        fprintf( stderr, "bad # args\n" );
        exit(1);
    }
    for ( int i=0; i<9; ++i )
        A[i] = atof( argv[1+i] );
    matrix_transposition( 3, A );

    printf( "Input matrix A:\n" );
    print_matrix(3, 3, A);

    lm_qrfac(3, 3, A, ipvt, rdiag, acnorm, wa);

    printf( "Output matrix A:\n" );
    print_matrix(3, 3, A);

    printf( "Output vector rdiag:\n" );
    print_matrix(1, 3, rdiag);

    printf( "Output vector acnorm:\n" );
    print_matrix(1, 3, acnorm);

    printf( "Output vector ipvt:\n" );
    for (i=0; i<3; ++i)
        printf( "%i ", ipvt[i] );
    printf("\n");
}