Eigen/src/SparseCholesky/SimplicialCholesky_impl.h - eigen/fuchsia - Git at Google

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr>

 /*

 NOTE: thes functions vave been adapted from the LDL library:

 LDL Copyright (c) 2005 by Timothy A. Davis.  All Rights Reserved.

 LDL License:

     Your use or distribution of LDL or any modified version of
     LDL implies that you agree to this License.

     This library is free software; you can redistribute it and/or
     modify it under the terms of the GNU Lesser General Public
     License as published by the Free Software Foundation; either
     version 2.1 of the License, or (at your option) any later version.

     This library is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     Lesser General Public License for more details.

     You should have received a copy of the GNU Lesser General Public
     License along with this library; if not, write to the Free Software
     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
     USA

     Permission is hereby granted to use or copy this program under the
     terms of the GNU LGPL, provided that the Copyright, this License,
     and the Availability of the original version is retained on all copies.
     User documentation of any code that uses this code or any modified
     version of this code must cite the Copyright, this License, the
     Availability note, and "Used by permission." Permission to modify
     the code and to distribute modified code is granted, provided the
     Copyright, this License, and the Availability note are retained,
     and a notice that the code was modified is included.
  */

 #include "../Core/util/NonMPL2.h"

 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
 #define EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H

 namespace Eigen {

 template<typename Derived>
 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
 {
   const StorageIndex size = StorageIndex(ap.rows());
   m_matrix.resize(size, size);
   m_parent.resize(size);
   m_nonZerosPerCol.resize(size);

   ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);

   for(StorageIndex k = 0; k < size; ++k)
   {
     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
     m_parent[k] = -1;             /* parent of k is not yet known */
     tags[k] = k;                  /* mark node k as visited */
     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
     for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
     {
       StorageIndex i = it.index();
       if(i < k)
       {
         /* follow path from i to root of etree, stop at flagged node */
         for(; tags[i] != k; i = m_parent[i])
         {
           /* find parent of i if not yet determined */
           if (m_parent[i] == -1)
             m_parent[i] = k;
           m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
           tags[i] = k;                  /* mark i as visited */
         }
       }
     }
   }

   /* construct Lp index array from m_nonZerosPerCol column counts */
   StorageIndex* Lp = m_matrix.outerIndexPtr();
   Lp[0] = 0;
   for(StorageIndex k = 0; k < size; ++k)
     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);

   m_matrix.resizeNonZeros(Lp[size]);

   m_isInitialized     = true;
   m_info              = Success;
   m_analysisIsOk      = true;
   m_factorizationIsOk = false;
 }


 template<typename Derived>
 template<bool DoLDLT>
 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
 {
   using std::sqrt;

   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
   eigen_assert(ap.rows()==ap.cols());
   eigen_assert(m_parent.size()==ap.rows());
   eigen_assert(m_nonZerosPerCol.size()==ap.rows());

   const StorageIndex size = StorageIndex(ap.rows());
   const StorageIndex* Lp = m_matrix.outerIndexPtr();
   StorageIndex* Li = m_matrix.innerIndexPtr();
   Scalar* Lx = m_matrix.valuePtr();

   ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
   ei_declare_aligned_stack_constructed_variable(StorageIndex,  pattern, size, 0);
   ei_declare_aligned_stack_constructed_variable(StorageIndex,  tags, size, 0);

   bool ok = true;
   m_diag.resize(DoLDLT ? size : 0);

   for(StorageIndex k = 0; k < size; ++k)
   {
     // compute nonzero pattern of kth row of L, in topological order
     y[k] = 0.0;                     // Y(0:k) is now all zero
     StorageIndex top = size;               // stack for pattern is empty
     tags[k] = k;                    // mark node k as visited
     m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
     for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
     {
       StorageIndex i = it.index();
       if(i <= k)
       {
         y[i] += numext::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
         Index len;
         for(len = 0; tags[i] != k; i = m_parent[i])
         {
           pattern[len++] = i;     /* L(k,i) is nonzero */
           tags[i] = k;            /* mark i as visited */
         }
         while(len > 0)
           pattern[--top] = pattern[--len];
       }
     }

     /* compute numerical values kth row of L (a sparse triangular solve) */

     RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset;    // get D(k,k), apply the shift function, and clear Y(k)
     y[k] = 0.0;
     for(; top < size; ++top)
     {
       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
       Scalar yi = y[i];             /* get and clear Y(i) */
       y[i] = 0.0;

       /* the nonzero entry L(k,i) */
       Scalar l_ki;
       if(DoLDLT)
         l_ki = yi / m_diag[i];
       else
         yi = l_ki = yi / Lx[Lp[i]];

       Index p2 = Lp[i] + m_nonZerosPerCol[i];
       Index p;
       for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
         y[Li[p]] -= numext::conj(Lx[p]) * yi;
       d -= numext::real(l_ki * numext::conj(yi));
       Li[p] = k;                          /* store L(k,i) in column form of L */
       Lx[p] = l_ki;
       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
     }
     if(DoLDLT)
     {
       m_diag[k] = d;
       if(d == RealScalar(0))
       {
         ok = false;                         /* failure, D(k,k) is zero */
         break;
       }
     }
     else
     {
       Index p = Lp[k] + m_nonZerosPerCol[k]++;
       Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
       if(d <= RealScalar(0)) {
         ok = false;              /* failure, matrix is not positive definite */
         break;
       }
       Lx[p] = sqrt(d) ;
     }
   }

   m_info = ok ? Success : NumericalIssue;
   m_factorizationIsOk = true;
 }

 } // end namespace Eigen

 #endif // EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra.
	//
	// Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr>

	/*

	NOTE: thes functions vave been adapted from the LDL library:

	LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.

	LDL License:

	Your use or distribution of LDL or any modified version of
	LDL implies that you agree to this License.

	This library is free software; you can redistribute it and/or
	modify it under the terms of the GNU Lesser General Public
	License as published by the Free Software Foundation; either
	version 2.1 of the License, or (at your option) any later version.

	This library is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public
	License along with this library; if not, write to the Free Software
	Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
	USA

	Permission is hereby granted to use or copy this program under the
	terms of the GNU LGPL, provided that the Copyright, this License,
	and the Availability of the original version is retained on all copies.
	User documentation of any code that uses this code or any modified
	version of this code must cite the Copyright, this License, the
	Availability note, and "Used by permission." Permission to modify
	the code and to distribute modified code is granted, provided the
	Copyright, this License, and the Availability note are retained,
	and a notice that the code was modified is included.
	*/

	#include "../Core/util/NonMPL2.h"

	#ifndef EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
	#define EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H

	namespace Eigen {

	template<typename Derived>
	void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
	{
	const StorageIndex size = StorageIndex(ap.rows());
	m_matrix.resize(size, size);
	m_parent.resize(size);
	m_nonZerosPerCol.resize(size);

	ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);

	for(StorageIndex k = 0; k < size; ++k)
	{
	/* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
	m_parent[k] = -1; /* parent of k is not yet known */
	tags[k] = k; /* mark node k as visited */
	m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
	for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
	{
	StorageIndex i = it.index();
	if(i < k)
	{
	/* follow path from i to root of etree, stop at flagged node */
	for(; tags[i] != k; i = m_parent[i])
	{
	/* find parent of i if not yet determined */
	if (m_parent[i] == -1)
	m_parent[i] = k;
	m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
	tags[i] = k; /* mark i as visited */
	}
	}
	}
	}

	/* construct Lp index array from m_nonZerosPerCol column counts */
	StorageIndex* Lp = m_matrix.outerIndexPtr();
	Lp[0] = 0;
	for(StorageIndex k = 0; k < size; ++k)
	Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);

	m_matrix.resizeNonZeros(Lp[size]);

	m_isInitialized = true;
	m_info = Success;
	m_analysisIsOk = true;
	m_factorizationIsOk = false;
	}


	template<typename Derived>
	template<bool DoLDLT>
	void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
	{
	using std::sqrt;

	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
	eigen_assert(ap.rows()==ap.cols());
	eigen_assert(m_parent.size()==ap.rows());
	eigen_assert(m_nonZerosPerCol.size()==ap.rows());

	const StorageIndex size = StorageIndex(ap.rows());
	const StorageIndex* Lp = m_matrix.outerIndexPtr();
	StorageIndex* Li = m_matrix.innerIndexPtr();
	Scalar* Lx = m_matrix.valuePtr();

	ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
	ei_declare_aligned_stack_constructed_variable(StorageIndex, pattern, size, 0);
	ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);

	bool ok = true;
	m_diag.resize(DoLDLT ? size : 0);

	for(StorageIndex k = 0; k < size; ++k)
	{
	// compute nonzero pattern of kth row of L, in topological order
	y[k] = 0.0; // Y(0:k) is now all zero
	StorageIndex top = size; // stack for pattern is empty
	tags[k] = k; // mark node k as visited
	m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
	for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
	{
	StorageIndex i = it.index();
	if(i <= k)
	{
	y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
	Index len;
	for(len = 0; tags[i] != k; i = m_parent[i])
	{
	pattern[len++] = i; /* L(k,i) is nonzero */
	tags[i] = k; /* mark i as visited */
	}
	while(len > 0)
	pattern[--top] = pattern[--len];
	}
	}

	/* compute numerical values kth row of L (a sparse triangular solve) */

	RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
	y[k] = 0.0;
	for(; top < size; ++top)
	{
	Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
	Scalar yi = y[i]; /* get and clear Y(i) */
	y[i] = 0.0;

	/* the nonzero entry L(k,i) */
	Scalar l_ki;
	if(DoLDLT)
	l_ki = yi / m_diag[i];
	else
	yi = l_ki = yi / Lx[Lp[i]];

	Index p2 = Lp[i] + m_nonZerosPerCol[i];
	Index p;
	for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
	y[Li[p]] -= numext::conj(Lx[p]) * yi;
	d -= numext::real(l_ki * numext::conj(yi));
	Li[p] = k; /* store L(k,i) in column form of L */
	Lx[p] = l_ki;
	++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
	}
	if(DoLDLT)
	{
	m_diag[k] = d;
	if(d == RealScalar(0))
	{
	ok = false; /* failure, D(k,k) is zero */
	break;
	}
	}
	else
	{
	Index p = Lp[k] + m_nonZerosPerCol[k]++;
	Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
	if(d <= RealScalar(0)) {
	ok = false; /* failure, matrix is not positive definite */
	break;
	}
	Lx[p] = sqrt(d) ;
	}
	}

	m_info = ok ? Success : NumericalIssue;
	m_factorizationIsOk = true;
	}

	} // end namespace Eigen

	#endif // EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H