testLapack.cpp 4.04 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================

#include <iostream>
#include <stdlib.h>

24 25
#include "Utils/FBlas.hpp"
#include "Utils/FTic.hpp"
26

27
#include "Utils/FParameterNames.hpp"
28

29 30 31 32
/**
 * Test functionality of C - interfaced LAPACK functions
 */

33
int main(int argc, char ** argv)
34
{
35
    FHelpDescribeAndExit(argc, argv, "Test the lapack compilation and linking (only the code is interesting).");
36

37
    /*
38 39 40 41 42 43
   * List of tested functions:
   * Cholesky decomposition: FBlas::potrf()
   * TODO SVD: FBlas::gesvd()
   * TODO QR decomposition: FBlas::geqrf()
   */

44
    typedef double FReal;
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
    const unsigned int m = 4, n = 4;
    FReal* A = new FReal [m * n]; // matrix: column major ordering

    // A= LL^T ////////////////////////////////////
    // define symmetric definite positive matrix A
    A[0]=5; A[10]=4; A[15]=7;
    A[1]=A[3]=A[4]=A[12]=2;
    A[6]=A[7]=A[9]=A[13]=1;
    A[2]=A[5]=A[8]=3;
    A[11]=A[14]=-1;

    // copy A in C
    FReal* C = new FReal [m * n]; // matrix: column major ordering
    for (unsigned int ii=0; ii<m; ++ii)
        for (unsigned int jj=0; jj<n; ++jj)
            C[ii*m + jj]=A[ii*m + jj];

    std::cout<<"\nA=["<<std::endl;
    for (unsigned int i=0; i<m; ++i) {
        for (unsigned int j=0; j<n; ++j)
            std::cout << A[i*n+j] << " ";
        std::cout<< std::endl;
67
    }
68 69 70 71 72 73 74 75 76 77 78 79 80
    std::cout<<"]"<<std::endl;

    // perform Cholesky decomposition
    std::cout<<"\nCholesky decomposition ";
    int INF = FBlas::potrf(m, A, n);
    if(INF==0) {std::cout<<"succeeded!"<<std::endl;}
    else {std::cout<<"failed!"<<std::endl;}

    std::cout<<"\nA_out=["<<std::endl;
    for (unsigned int i=0; i<m; ++i) {
        for (unsigned int j=0; j<n; ++j)
            std::cout << A[i*n+j] << " ";
        std::cout<<std::endl;
81
    }
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
    std::cout<<"]"<<std::endl;

    // build lower matrix
    FReal* L = new FReal [m * n]; // matrix: column major ordering
    for (unsigned int ii=0; ii<m; ++ii)
        for (unsigned int jj=0; jj<n; ++jj){
            if(ii<=jj)
                L[ii*m + jj]=A[ii*m + jj];
            else
                L[ii*m + jj]=0.;
        }

    std::cout<<"\nL=["<<std::endl;
    for (unsigned int i=0; i<m; ++i) {
        for (unsigned int j=0; j<n; ++j)
            std::cout << L[i*n+j] << " ";
        std::cout<< std::endl;
    }
    std::cout<<"]"<<std::endl;

    // verify result by computing B=LL^T
    FReal* B = new FReal [m * n]; // matrix: column major ordering
    for (unsigned int ii=0; ii<m; ++ii)
        for (unsigned int jj=0; jj<n; ++jj){
            B[ii*m + jj]=0.;
            for (unsigned int j=0; j<n; ++j)
                B[ii*m + jj]+=L[j*m + ii]*L[j*m + jj];
        }

    std::cout<<"\nA-LL^T=["<<std::endl;
    for (unsigned int i=0; i<m; ++i) {
        for (unsigned int j=0; j<n; ++j)
            std::cout << B[i*n+j]-C[i*n+j] << " ";
        std::cout<< std::endl;
    }
    std::cout<<"]"<<std::endl;
118

119 120 121
    delete [] A;
    delete [] B;
    delete [] C;
122
    delete [] L;
123

124
    return 0;
125
}