Commit 206ffcac authored by Peter Bastian's avatar Peter Bastian
Browse files

starte qr und gram schmidt

parent bc241986

examples/num0/integralgleichung.cc

+17 −0
Original line number Diff line number Diff line
@@ -82,6 +82,23 @@ int main () 
  // Diskretisierungsparameter

  int n = 10;



  // Teste Gram-Schmidt Orthogonalisierung

  {

    n = 128;

    gamma = 1e-4;

    DenseMatrix<double> A_double(setupA(integrateK,alpha,beta,n));

    Vector<double> b_double(setupb(f,alpha,beta,n));

    DenseMatrix<double> Q1 = gram_schmidt(A_double);

    DenseMatrix<double> Q1T(Q1.transpose());

    DenseMatrix<double> I1(A_double); I1.mm(Q1T,Q1);

    for (int i=0; i<Q1T.rowsize(); i++) I1[i][i] -= 1.0;

    DenseMatrix<double> Q2 = modified_gram_schmidt(A_double);

    DenseMatrix<double> Q2T(Q2.transpose());

    DenseMatrix<double> I2(A_double); I2.mm(Q2T,Q2);

    for (int i=0; i<Q2T.rowsize(); i++) I2[i][i] -= 1.0;

    std::cout << "cgs " << I1.norm_infty() << " mgs " << I2.norm_infty() << std::endl;

  }

  

  std::vector<int> sizes = {8,32,128,512};

  std::vector<double> gammas = {0.9,0.5,1e-1,1e-2,1e-3,1e-4,1e-5};

hdnum.hh

+1 −0
Original line number Diff line number Diff line
@@ -23,6 +23,7 @@ 
#include "src/sgrid.hh"

#include "src/newton.hh"

#include "src/rungekutta.hh"

#include "src/qr.hh"





#endif

src/qr.hh

0 → 100644
+91 −0
Original line number Diff line number Diff line 
// -*- tab-width: 4; indent-tabs-mode: nil -*-

#ifndef HDNUM_QR_HH

#define HDNUM_QR_HH



#include "vector.hh"

#include "densematrix.hh"

#include <cmath>



/** @file

 *  @brief This file implements QR decomposition

 */



namespace hdnum

{

  //! computes orthonormal basis of Im(A) using classical Gram-Schmidt

  template<class T>

  DenseMatrix<T> gram_schmidt (const DenseMatrix<T>& A)

  {

    DenseMatrix<T> Q(A);



    // for all columns except the first

    for (int k=1; k<Q.colsize(); k++)

      {

        // orthogonalize column k against all previous

        for (int j=0; j<k; j++)

          {

            // compute factor

            T sum_nom(0.0);

            T sum_denom(0.0);

            for (int i=0; i<Q.rowsize(); i++)

              {

                sum_nom += A[i][k]*Q[i][j];

                sum_denom += Q[i][j]*Q[i][j];

              }

            // modify 

            T alpha = sum_nom/sum_denom;

            for (int i=0; i<Q.rowsize(); i++)

              Q[i][k] -= alpha*Q[i][j];

          }

      }

    for (int j=0; j<Q.colsize(); j++)

      {

        // compute norm of column j

        T sum(0.0);

        for (int i=0; i<Q.rowsize(); i++) sum += Q[i][j]*Q[i][j];

        sum = sqrt(sum);

        //scale

        for (int i=0; i<Q.rowsize(); i++) Q[i][j] = Q[i][j]/sum;

      }

    return Q;

  }



  //! computes orthonormal basis of Im(A) using modified Gram-Schmidt 

  template<class T>

  DenseMatrix<T> modified_gram_schmidt (const DenseMatrix<T>& A)

  {

    DenseMatrix<T> Q(A);



    for (int k=0; k<Q.colsize(); k++)

      {

        // modify all later columns with column k

        for (int j=k+1; j<Q.rowsize(); j++)

          {

            // compute factor

            T sum_nom(0.0);

            T sum_denom(0.0);

            for (int i=0; i<Q.rowsize(); i++)

              {

                sum_nom += Q[i][j]*Q[i][k];

                sum_denom += Q[i][k]*Q[i][k];

              }

            // modify 

            T alpha = sum_nom/sum_denom;

            for (int i=0; i<Q.rowsize(); i++)

              Q[i][j] -= alpha*Q[i][k]; 

          }

      }

    for (int j=0; j<Q.colsize(); j++)

      {

        // compute norm of column j

        T sum(0.0);

        for (int i=0; i<Q.rowsize(); i++) sum += Q[i][j]*Q[i][j];

        sum = sqrt(sum);

        //scale

        for (int i=0; i<Q.rowsize(); i++) Q[i][j] = Q[i][j]/sum;

      }    

    return Q;

  }



}

#endif