mattest.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
 
#include <sp/spLib.h>
#include <sp/spMain.h>
 
int spMain(int argc, char **argv)
{
    long i;
    double det, detr, deti;
    spDMatrix mat, mat2;
    spDMatrix matinv;
    spDMatrix eye;
    spLVector index;
    spDMatrix L, U, D, V;
    spDMatrix B;
    spDMatrix Q, R, QR;
    spDVector b, x, y, v;
    spDVectors eigvecs;
 
    if (0 && argc >= 2) {
        spDMatrix S;
        
        mat = xdmreaddmatrix(argv[1], 6, 0);
        //dmdump(mat);
        S = xdmsvd(mat, 0, 0.0, NULL, NULL);
        dmdump(S);
        xdmfree(S);
        printf("cond = %g\n", dmcond(mat));
 
        S = xdminv(mat);
        dmdump(S);
        xdmfree(S);
        
        xdmfree(mat);
        spExit(0);
    }
    
    printf("---- Matrix No.1 ----\n");
#if 1
    mat = xdmzeros(3, 3);
#if 0
    mat->data[0][1] = 1;
    mat->data[1][2] = 1;
    mat->data[2][0] = 1;
#else
    mat->data[0][0] = 1;
    mat->data[0][1] = 2;
    mat->data[0][2] = 3;
#if 0
    mat->data[1][0] = 4;
    mat->data[1][1] = 5;
    mat->data[1][2] = 6;
    mat->data[2][0] = 7;
    mat->data[2][1] = 8;
    mat->data[2][2] = 0;
#elif 0
    mat->data[1][0] = 7;
    mat->data[1][1] = 8;
    mat->data[1][2] = 0;
    mat->data[2][0] = 4;
    mat->data[2][1] = 5;
    mat->data[2][2] = 6;
#else
    /* complex eigenvalues */
    mat->data[1][0] = 3;
    mat->data[1][1] = 1;
    mat->data[1][2] = 2;
    mat->data[2][0] = 2;
    mat->data[2][1] = 3;
    mat->data[2][2] = 1;
#endif
#endif
#elif 0
    mat = xdmzeros(4, 4);
#if 1
    /* complex eigenvalues */
    mat->data[0][1] = 1;
    mat->data[1][2] = 1;
    mat->data[2][3] = 1;
    mat->data[3][0] = 1;
#else
    /* complex eigenvalues; matlab rsf2csf sample */
    mat->data[0][0] = 1;
    mat->data[0][1] = 1;
    mat->data[0][2] = 1;
    mat->data[0][3] = 3;
    mat->data[1][0] = 1;
    mat->data[1][1] = 2;
    mat->data[1][2] = 1;
    mat->data[1][3] = 1;
    mat->data[2][0] = 1;
    mat->data[2][1] = 1;
    mat->data[2][2] = 3;
    mat->data[2][3] = 1;
    mat->data[3][0] = -2;
    mat->data[3][1] = 1;
    mat->data[3][2] = 1;
    mat->data[3][3] = 4;
#endif
#else
    mat = xdmzeros(4, 2);
    mat->data[0][0] = 1;
    mat->data[0][1] = 2;
    mat->data[1][0] = 3;
    mat->data[1][1] = 4;
    mat->data[2][0] = 5;
    mat->data[2][1] = 6;
    mat->data[3][0] = 7;
    mat->data[3][1] = 8;
#endif
 
    eye = xdmeye(mat->row, mat->col);
 
    dmdump(mat);
    dmdump(eye);
 
    printf("cond of mat = %f\n", dmcond(mat));
    printf("cond of eye = %f\n", dmcond(eye));
 
    /* QR decomposition */
    Q = xdmalloc(mat->row, mat->row);
    R = xdmalloc(mat->row, mat->col);
 
    dmqr(mat, Q, R);
    printf("Q of QR decomposition:\n");
    dmdump(Q);
    printf("R of QR decomposition:\n");
    dmdump(R);
 
    QR = xdmtimes(Q, R);
    printf("Q*R:\n");
    dmdump(QR);
 
    xdmfree(Q);
    xdmfree(R);
    xdmfree(QR);
 
    if ((B = xdmbidiag(mat, &U, &V)) != NODATA) {
        spDMatrix UB, UBV;
        printf("bidiag result:\n");
        dmdump(B);
        dmdump(U);
        dmdump(V);
 
        UB = xdmtimes(U, B);
        dmconjtranspose(V);
        UBV = xdmtimes(UB, V);
        dmdump(UBV);
        
        xdmfree(B);
        xdmfree(U);
        xdmfree(V);
        xdmfree(UB);
        xdmfree(UBV);
    }
 
    if ((B = xdmsvd(mat, 50, 0.0, &U, &V)) != NODATA) {
        spDMatrix UB, UBV;
        printf("SVD result:\n");
        dmdump(B);
        dmdump(U);
        dmdump(V);
 
        UB = xdmtimes(U, B);
        dmconjtranspose(V);
        UBV = xdmtimes(UB, V);
        dmdump(UBV);
        
        xdmfree(B);
        xdmfree(U);
        xdmfree(V);
        xdmfree(UB);
        xdmfree(UBV);
    }
 
    if ((y = xdmeigdsqr(mat, 50, 0.0, &eigvecs)) != NODATA) {
        printf("eigdsqr result:\n");
        dvdump(y);
        for (i = 0; i < y->length; i++) {
            if (y->imag != NULL) {
                printf("eigenvector for %f + %fi:\n", y->data[i], y->imag[i]);
            } else {
                printf("eigenvector for %f:\n", y->data[i]);
            }
            if (eigvecs->vector[i] != NODATA) {
                dvdump(eigvecs->vector[i]);
            } else {
                printf("Can't get eigenvector: No. %ld\n", i);
            }
        }
        xdvsfree(eigvecs);
        xdvfree(y);
    }
 
    /* LU decomposition */
    if (dmlu(mat, &index, &det) == SP_TRUE) {
        printf("det = %f\n", det);
 
        L = xdmlulower(mat, index);
        U = xdmluupper(mat, index);
        D = xdmludiag(mat, index);
 
        matinv = xdmluinv(mat, index);
 
        dmdump(mat);
        lvdump(index);
        /*dmdump(eye);*/
        dmdump(L);
        dmdump(U);
        dmdump(D);
        dmdump(matinv);
    
        xdmfree(L);
        xdmfree(U);
        xdmfree(D);
        xdmfree(eye);
        xdmfree(matinv);
    
        b = xdvinit(4.0, 2.0, 8.0);
        y = xdvlusolve(mat, index, b);
        dvdump(b);
        dvdump(y);
 
        xlvfree(index);
        xdvfree(y);
    
        xdvfree(b);
    }
 
    xdmfree(mat);
 
    mat = xdmpascal(5);
    b = xdvinit(4.0, 2.0, 12.0);
    L = xdmchol(mat, SP_TRUE);
    U = xdmchol(mat, SP_FALSE);
    dmdump(mat);
    
    dmdump(L);
    dmdump(U);
    
    mat2 = xdmtimes(L, U);
    dmdump(mat2);
    xdmfree(mat2);
 
    printf("cholsolve (lower) result:\n");
    y = xdvcholsolve(L, b, SP_TRUE);
    dvdump(y);
    xdvfree(y);
    
    printf("cholsolve (uppwer) result:\n");
    y = xdvcholsolve(U, b, SP_FALSE);
    dvdump(y);
    xdvfree(y);
    
    xdmfree(L);
    xdmfree(U);
 
    mat->data[4][4] -= 1.0;     /* make non-positive definite */
    if ((L = xdmchol(mat, SP_TRUE)) == NODATA) {
        printf("xdmchol error: matrix must be positive definite:\n");
    } else {
        dmdump(L);
        printf("xdmchol problem: non-positive definite matrix has been passed:\n");
        xdmfree(L);
    }
    dmdump(mat);
    xdmfree(mat);
 
    printf("---- Matrix No.2 ----\n");
    mat = xdmrizeros(3, 3);
#if 1
    dmiones(mat, 0, 0);
    mat->data[0][0] = -2.0;
    mat->data[1][1] = 2.0;
    mat->imag[1][1] = 0.0;
    mat->data[2][0] = 1.0;
    mat->data[2][2] = 2.0;
#if 0
    mat->imag[2][2] = 0.0;
#else
    mat->imag[2][2] = 0.5;
#endif
#elif 1
    /* has complex conjugate eigenvalues */
    mat->data[0][1] = 1.0;
    mat->data[1][2] = 1.0;
    mat->data[2][0] = 1.0;
#else
    mat->data[0][0] = 1.0;
    mat->data[0][1] = 100.0;
    mat->data[0][2] = 10000.0;
    mat->data[1][0] = 0.01;
    mat->data[1][1] = 1.0;
    mat->data[1][2] = 100.0;
    mat->data[2][0] = 0.0001;
    mat->data[2][1] = 0.01;
    mat->data[2][2] = 1.0;
#endif
    
    dmdump(mat);
 
    /* QR decomposition */
    Q = xdmalloc(mat->row, mat->row);
    R = xdmalloc(mat->row, mat->col);
 
    dmqr(mat, Q, R);
    printf("Q of QR decomposition:\n");
    dmdump(Q);
    printf("R of QR decomposition:\n");
    dmdump(R);
 
    QR = xdmtimes(Q, R);
    printf("Q*R:\n");
    dmdump(QR);
 
    xdmfree(Q);
    xdmfree(R);
    xdmfree(QR);
 
    if ((B = xdmbidiag(mat, &U, &V)) != NODATA) {
        spDMatrix UB, UBV;
        printf("bidiag result:\n");
        dmdump(B);
        dmdump(U);
        dmdump(V);
 
        UB = xdmtimes(U, B);
        dmconjtranspose(V);
        UBV = xdmtimes(UB, V);
        dmdump(UBV);
        
        xdmfree(B);
        xdmfree(U);
        xdmfree(V);
        xdmfree(UB);
        xdmfree(UBV);
    }
    
    if ((B = xdmsvd(mat, 50, 0.0, &U, &V)) != NODATA) {
        spDMatrix UB, UBV;
        printf("SVD result:\n");
        dmdump(B);
        dmdump(U);
        dmdump(V);
 
        UB = xdmtimes(U, B);
        dmconjtranspose(V);
        UBV = xdmtimes(UB, V);
        dmdump(UBV);
        
        xdmfree(B);
        xdmfree(U);
        xdmfree(V);
        xdmfree(UB);
        xdmfree(UBV);
    }
    
    /* calculate eigenvector */
    if (0) {
        dmbalance(mat, NODATA);
        dmdump(mat);
    }
    if ((y = xdmeigdsqr(mat, 50, 0.0, &eigvecs)) != NODATA) {
        printf("eigdsqr result:\n");
        dvdump(y);
        for (i = 0; i < y->length; i++) {
            if (y->imag != NULL) {
                printf("eigenvector for %f + %fi:\n", y->data[i], y->imag[i]);
            } else {
                printf("eigenvector for %f:\n", y->data[i]);
            }
            if (eigvecs->vector[i] != NODATA) {
                dvdump(eigvecs->vector[i]);
            } else {
                printf("Can't get eigenvector: No. %ld\n", i);
            }
        }
        xdvfree(y);
    }
 
    /* Cholesky decomposition for non-Hermitian (wrong case just for check) */
    L = xdmchol(mat, SP_TRUE);
    U = xdmchol(mat, SP_FALSE);
    if (L != NODATA && U != NODATA) {
        mat2 = xdmtimes(L, U);
        dmdump(mat2);
        xdmfree(mat2);
    
        dmdump(L);
        dmdump(U);
    }
 
    if (L != NODATA) {
        printf("cholinv (lower) result:\n");
        matinv = xdmcholinv(L, SP_TRUE);
        dmdump(matinv);
        xdmfree(matinv);
    }
    
    if (U != NODATA) {
        printf("cholinv (upper) result:\n");
        matinv = xdmcholinv(U, SP_FALSE);
        dmdump(matinv);
        xdmfree(matinv);
    }
    
    if (dmlucplx(mat, &index, &detr, &deti) == SP_TRUE) {
        printf("det = %f + %fj\n", detr, deti);
        dmdump(mat);
        matinv = xdmluinv(mat, index);
        dmdump(matinv);
        b = xdvinit(4.0, 2.0, 8.0);
        y = xdvlusolve(mat, index, b);
        printf("lusolve result:\n");
        dvdump(y);
        xdvfree(y);
        y = xdvmvtimes(matinv, b);
        dvdump(y);
        xdmfree(matinv);
        xdvfree(y);
        //lvdump(index);
 
        if (L != NODATA) {
            printf("cholsolve (lower) result:\n");
            y = xdvcholsolve(L, b, SP_TRUE);
            dvdump(y);
            xdvfree(y);
        }
 
        if (U != NODATA) {
            printf("cholsolve (uppwer) result:\n");
            y = xdvcholsolve(U, b, SP_FALSE);
            dvdump(y);
            xdvfree(y);
        }
    
        xdvfree(b);
        xlvfree(index);
    }
    
    xdmfree(L);
    xdmfree(U);
    
    xdmfree(mat);
 
    /* householder matrix */
    x = xdvinit(1, 1, 4);
    printf("householder matrix test input vector:\n");
    dvdump(x);
    v = xdvhouse(x, NULL);
    printf("householder hyperplane normal vector:\n");
    dvdump(v);
    printf("householder matrix:\n");
    mat = xdmhouse(v);
    dmdump(mat);
    y = xdvmvtimes(mat, x);
    printf("converted vector by householder matrix:\n");
    dvdump(y);
    xdvfree(x);
    xdvfree(v);
    xdvfree(y);
    xdmfree(mat);
 
    x = xdvriinit(1, 1, 4);
    //x->imag[0] *= -1.0;
    x->imag[2] *= -1.0;
    //x = xdvinit(1, 1, 3);
    printf("householder matrix test input vector (complex case):\n");
    dvdump(x);
    v = xdvhouse(x, NULL);
    printf("householder hyperplane normal vector (complex case):\n");
    dvdump(v);
    mat = xdmhouse(v);
    dmdump(mat);
    y = xdvmvtimes(mat, x);
    printf("converted vector by householder matrix (complex case):\n");
    dvdump(y);
    xdvfree(x);
    xdvfree(v);
    xdvfree(y);
    xdmfree(mat);
 
#if 0
    v = xdvzeros(3);
    v->data[0] = 3.0;
    v->data[1] = -2.0;
    v->data[2] = -4.0;
#else
    v = xdvzeros(4);
    v->data[0] = 1.0;
    v->data[2] = -7.0;
    v->data[3] = 6.0;
#endif
    printf("compute roots for:\n");
    dvdump(v);
    if ((mat = xdmcompan(v, 0.0, NULL, NULL)) != NODATA) {
        if ((x = xdvrootscompan(mat, 0, 0)) != NODATA) {
            printf("roots result:\n");
            dvdump(x);
            xdvfree(x);
        }
        printf("companion matrix:\n");
        dmdump(mat);
        
        xdmfree(mat);
    }
    xdvfree(v);
    
    v = xdvzeros(5);
    v->data[0] = 1.0;
    v->data[4] = -1.0;
    printf("compute roots for:\n");
    dvdump(v);
    if ((mat = xdmcompan(v, 0.0, NULL, NULL)) != NODATA) {
        if ((x = xdvrootscompan(mat, 0, 0)) != NODATA) {
            printf("roots result:\n");
            dvdump(x);
            xdvfree(x);
        }
        printf("companion matrix:\n");
        dmdump(mat);
        
        xdmfree(mat);
    }
    xdvfree(v);
    
    return 0;
}