testLapackQR.cpp 4.57 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
// ===================================================================================
// 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>
BLANCHARD Pierre's avatar
BLANCHARD Pierre committed
23
#include <stdexcept>
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 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 118

#include "Utils/FGlobal.hpp"
#include "Utils/FBlas.hpp"
#include "Utils/FTic.hpp"
#include "Utils/FParameterNames.hpp"

/**
 * Test functionality of C - interfaced LAPACK functions : GEQP3
 */

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

    /*
   * List of tested functions:
   * QR decomposition with pivoting: FBlas::geqp3()
   */
    
    FTic time;

    typedef double FReal;
    const FReal FRandMax = FReal(RAND_MAX);

    const unsigned int m = 10, n = 10;
    FReal* A = new FReal [m * n]; // matrix: column major ordering

    // A= LL^T ////////////////////////////////////
    // define symmetric definite positive matrix A
    for (unsigned int i=0; i<m; ++i) 
        for (unsigned int j=0; j<n; ++j)
            A[i+j*m]=FReal(rand())/FRandMax;

    // 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 + jj*m]=A[ii + jj*m];

    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+j*m] << " ";
        std::cout<< std::endl;
    }
    std::cout<<"]"<<std::endl;

    // Workspace
    // init SVD
    const unsigned int minMN = std::min(m,n);
    const unsigned int maxMN = std::max(m,n);
    //const FSize LWORK = 2*4*minMN; // for square matrices
    const unsigned int LWORK = 2*std::max(3*minMN+maxMN, 5*minMN);
    FReal *const WORK = new FReal[LWORK];
    FReal* TAU = new FReal[n];

    // Pivot
    unsigned* jpiv = new unsigned[n]; 

    // perform Cholesky decomposition
    std::cout<<"\nQR decomposition ";
    FTic timeQR;
    //int INF = FBlas::geqrf(m, n, A, TAU, LWORK, WORK);
    int INFO = FBlas::geqp3(m, n, A, jpiv, TAU, LWORK, WORK);
     //FBlas::geqp3(m, n, A, jpiv);

    if(INFO==0) {std::cout<<"succeeded!"<<std::endl;}
    else {std::cout<<"failed!"<<std::endl;}



    std::cout << "took " << timeQR.tacAndElapsed() << "sec." << 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+j*m] << " ";
        std::cout<<std::endl;
    }
    std::cout<<"]"<<std::endl;

    std::cout<<"\njpiv=["<<std::endl;
    for (unsigned int i=0; i<n; ++i) {
            std::cout << jpiv[i] << " ";
    }
    std::cout<<"]"<<std::endl;


    // get Q
    FTic timeGETQ;
    timeGETQ.tic();
    INFO = FBlas::orgqr(m, n, A, TAU, int(LWORK), WORK);
    if(INFO!=0) {
      std::stringstream stream;
      stream << INFO;
119 120 121 122 123
      delete [] A;
      delete [] C;
      delete [] jpiv;
      delete [] TAU;
      delete [] WORK;
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
      throw std::runtime_error("get Q failed! INFO=" + stream.str());
    }
    double tGETQ = timeGETQ.tacAndElapsed();
    std::cout << "... took @tGETQ = "<< tGETQ <<"\n";


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

    delete [] TAU;
    delete [] WORK;

    // TODO Verify QR facto! How? Form R then apply Q? How?

    delete [] A;
    delete [] C;
    delete [] jpiv;

    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    return 0;
}