~ubuntu-branches/ubuntu/saucy/blender/saucy-proposed

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
4
4
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5
5
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6
6
//
7
 
// Eigen is free software; you can redistribute it and/or
8
 
// modify it under the terms of the GNU Lesser General Public
9
 
// License as published by the Free Software Foundation; either
10
 
// version 3 of the License, or (at your option) any later version.
11
 
//
12
 
// Alternatively, you can redistribute it and/or
13
 
// modify it under the terms of the GNU General Public License as
14
 
// published by the Free Software Foundation; either version 2 of
15
 
// the License, or (at your option) any later version.
16
 
//
17
 
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18
 
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19
 
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20
 
// GNU General Public License for more details.
21
 
//
22
 
// You should have received a copy of the GNU Lesser General Public
23
 
// License and a copy of the GNU General Public License along with
24
 
// Eigen. If not, see <http://www.gnu.org/licenses/>.
 
7
// This Source Code Form is subject to the terms of the Mozilla
 
8
// Public License v. 2.0. If a copy of the MPL was not distributed
 
9
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
25
10
 
26
11
#ifndef EIGEN_REAL_SCHUR_H
27
12
#define EIGEN_REAL_SCHUR_H
28
13
 
29
 
#include "./EigenvaluesCommon.h"
30
14
#include "./HessenbergDecomposition.h"
31
15
 
 
16
namespace Eigen { 
 
17
 
32
18
/** \eigenvalues_module \ingroup Eigenvalues_Module
33
19
  *
34
20
  *
235
221
  // Rows iu+1,...,end are already brought in triangular form.
236
222
  Index iu = m_matT.cols() - 1;
237
223
  Index iter = 0; // iteration count
238
 
  Scalar exshift = 0.0; // sum of exceptional shifts
 
224
  Scalar exshift(0); // sum of exceptional shifts
239
225
  Scalar norm = computeNormOfT();
240
226
 
241
 
  while (iu >= 0)
 
227
  if(norm!=0)
242
228
  {
243
 
    Index il = findSmallSubdiagEntry(iu, norm);
244
 
 
245
 
    // Check for convergence
246
 
    if (il == iu) // One root found
247
 
    {
248
 
      m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
249
 
      if (iu > 0) 
250
 
        m_matT.coeffRef(iu, iu-1) = Scalar(0);
251
 
      iu--;
252
 
      iter = 0;
253
 
    }
254
 
    else if (il == iu-1) // Two roots found
255
 
    {
256
 
      splitOffTwoRows(iu, computeU, exshift);
257
 
      iu -= 2;
258
 
      iter = 0;
259
 
    }
260
 
    else // No convergence yet
261
 
    {
262
 
      // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
263
 
      Vector3s firstHouseholderVector(0,0,0), shiftInfo;
264
 
      computeShift(iu, iter, exshift, shiftInfo);
265
 
      iter = iter + 1; 
266
 
      if (iter > m_maxIterations) break;
267
 
      Index im;
268
 
      initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
269
 
      performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
270
 
    }
271
 
  } 
272
 
 
273
 
  if(iter <= m_maxIterations) 
 
229
    while (iu >= 0)
 
230
    {
 
231
      Index il = findSmallSubdiagEntry(iu, norm);
 
232
 
 
233
      // Check for convergence
 
234
      if (il == iu) // One root found
 
235
      {
 
236
        m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
 
237
        if (iu > 0)
 
238
          m_matT.coeffRef(iu, iu-1) = Scalar(0);
 
239
        iu--;
 
240
        iter = 0;
 
241
      }
 
242
      else if (il == iu-1) // Two roots found
 
243
      {
 
244
        splitOffTwoRows(iu, computeU, exshift);
 
245
        iu -= 2;
 
246
        iter = 0;
 
247
      }
 
248
      else // No convergence yet
 
249
      {
 
250
        // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
 
251
        Vector3s firstHouseholderVector(0,0,0), shiftInfo;
 
252
        computeShift(iu, iter, exshift, shiftInfo);
 
253
        iter = iter + 1;
 
254
        if (iter > m_maxIterations * m_matT.cols()) break;
 
255
        Index im;
 
256
        initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
 
257
        performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
 
258
      }
 
259
    }
 
260
  }
 
261
  if(iter <= m_maxIterations * m_matT.cols()) 
274
262
    m_info = Success;
275
263
  else
276
264
    m_info = NoConvergence;
288
276
  // FIXME to be efficient the following would requires a triangular reduxion code
289
277
  // Scalar norm = m_matT.upper().cwiseAbs().sum() 
290
278
  //               + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
291
 
  Scalar norm = 0.0;
 
279
  Scalar norm(0);
292
280
  for (Index j = 0; j < size; ++j)
293
281
    norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
294
282
  return norm;
471
459
  }
472
460
}
473
461
 
 
462
} // end namespace Eigen
 
463
 
474
464
#endif // EIGEN_REAL_SCHUR_H