BFGS_Precond.cc 1.43 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
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
#include "BFGS_Precond.h"
#include "QN_Precond.h"
#include "FiniteElemSpace.h"
#include "DOFAdmin.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
#include "DOFIterator.h"

namespace AMDiS
{

  void BFGS_Precond::addPair(DOFVector<double> *s, DOFVector<double> *y,
			     double sy)
  {
    FUNCNAME("BFGS_Preconditioner::addPair()");

    if (!sy)
      sy = *s * *y;
    TEST_EXIT(sy > 0)("Vectors do not satisfy the curvature condition!\n");

    BFGS_Vectors newPair(s->getFESpace());
    newPair.setVectors(s, y, sy);

    if (k < m) {
      k++;
    } else {
      pairsVec.pop_back();
    }

    pairsVec.push_front(newPair);

    return;
  }

  void BFGS_Precond::precon(DOFVector<double> *vec)
  {
    FUNCNAME("BFGS_Preconditioner::precon()");

    int    i      = 0;
    double beta;
    double *alpha = GET_MEMORY(double, k);

    ::std::deque<BFGS_Vectors>::iterator pairsIterator(pairsVec.begin());

    while (pairsIterator != pairsVec.end()) {
      alpha[i]  = pairsIterator->s * *vec;
      alpha[i] /= pairsIterator->sy;
      aXpY(-alpha[i], pairsIterator->y, *vec, bound);
      ++i;
      ++pairsIterator;
    }

    diagPrecond(*matrix[row], vec, bound);

    while (pairsIterator != pairsVec.begin()) {
      --pairsIterator;
      --i;
      beta  = pairsIterator->y * *vec;
      beta /= pairsIterator->sy;
      aXpY(alpha[i]-beta, pairsIterator->s, *vec, bound);
    }
    
    FREE_MEMORY(alpha, double, k);
    
    return;
  }

}