src/parsec/disk-image/parsec/parsec-benchmark/pkgs/libs/gsl/src/linalg/cholesky.c - public/gem5-resources - Git at Google

 /* Cholesky Decomposition
  *
  * Copyright (C) 2000  Thomas Walter
  *
  * 03 May 2000: Modified for GSL by Brian Gough
  * 29 Jul 2005: Additions by Gerard Jungman
  * Copyright (C) 2000,2001, 2002, 2003, 2005 Brian Gough, Gerard Jungman
  *
  * This is free software; you can redistribute it and/or modify it
  * under the terms of the GNU General Public License as published by the
  * Free Software Foundation; either version 2, or (at your option) any
  * later version.
  *
  * This source 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 General Public License
  * for more details.
  */

 /*
  * Cholesky decomposition of a symmetrix positive definite matrix.
  * This is useful to solve the matrix arising in
  *    periodic cubic splines
  *    approximating splines
  *
  * This algorithm does:
  *   A = L * L'
  * with
  *   L  := lower left triangle matrix
  *   L' := the transposed form of L.
  *
  */

 #include <config.h>

 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_linalg.h>

 static inline
 double
 quiet_sqrt (double x)
      /* avoids runtime error, for checking matrix for positive definiteness */
 {
   return (x >= 0) ? sqrt(x) : GSL_NAN;
 }

 int
 gsl_linalg_cholesky_decomp (gsl_matrix * A)
 {
   const size_t M = A->size1;
   const size_t N = A->size2;

   if (M != N)
     {
       GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
     }
   else
     {
       size_t i,j,k;
       int status = 0;

       /* Do the first 2 rows explicitly.  It is simple, and faster.  And
        * one can return if the matrix has only 1 or 2 rows.
        */

       double A_00 = gsl_matrix_get (A, 0, 0);

       double L_00 = quiet_sqrt(A_00);

       if (A_00 <= 0)
         {
           status = GSL_EDOM ;
         }

       gsl_matrix_set (A, 0, 0, L_00);

       if (M > 1)
         {
           double A_10 = gsl_matrix_get (A, 1, 0);
           double A_11 = gsl_matrix_get (A, 1, 1);

           double L_10 = A_10 / L_00;
           double diag = A_11 - L_10 * L_10;
           double L_11 = quiet_sqrt(diag);

           if (diag <= 0)
             {
               status = GSL_EDOM;
             }

           gsl_matrix_set (A, 1, 0, L_10);
           gsl_matrix_set (A, 1, 1, L_11);
         }

       for (k = 2; k < M; k++)
         {
           double A_kk = gsl_matrix_get (A, k, k);

           for (i = 0; i < k; i++)
             {
               double sum = 0;

               double A_ki = gsl_matrix_get (A, k, i);
               double A_ii = gsl_matrix_get (A, i, i);

               gsl_vector_view ci = gsl_matrix_row (A, i);
               gsl_vector_view ck = gsl_matrix_row (A, k);

               if (i > 0) {
                 gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
                 gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);

                 gsl_blas_ddot (&di.vector, &dk.vector, &sum);
               }

               A_ki = (A_ki - sum) / A_ii;
               gsl_matrix_set (A, k, i, A_ki);
             }

           {
             gsl_vector_view ck = gsl_matrix_row (A, k);
             gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);

             double sum = gsl_blas_dnrm2 (&dk.vector);
             double diag = A_kk - sum * sum;

             double L_kk = quiet_sqrt(diag);

             if (diag <= 0)
               {
                 status = GSL_EDOM;
               }

             gsl_matrix_set (A, k, k, L_kk);
           }
         }

       /* Now copy the transposed lower triangle to the upper triangle,
        * the diagonal is common.
        */

       for (i = 1; i < M; i++)
         {
           for (j = 0; j < i; j++)
             {
               double A_ij = gsl_matrix_get (A, i, j);
               gsl_matrix_set (A, j, i, A_ij);
             }
         }

       if (status == GSL_EDOM)
         {
           GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
         }

       return GSL_SUCCESS;
     }
 }


 int
 gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
                            const gsl_vector * b,
                            gsl_vector * x)
 {
   if (LLT->size1 != LLT->size2)
     {
       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
     }
   else if (LLT->size1 != b->size)
     {
       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
     }
   else if (LLT->size2 != x->size)
     {
       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
     }
   else
     {
       /* Copy x <- b */

       gsl_vector_memcpy (x, b);

       /* Solve for c using forward-substitution, L c = b */

       gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

       /* Perform back-substitution, U x = c */

       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);


       return GSL_SUCCESS;
     }
 }

 int
 gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
                          gsl_vector * x)
 {
   if (LLT->size1 != LLT->size2)
     {
       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
     }
   else if (LLT->size2 != x->size)
     {
       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
     }
   else
     {
       /* Solve for c using forward-substitution, L c = b */

       gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

       /* Perform back-substitution, U x = c */

       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

       return GSL_SUCCESS;
     }
 }


 int
 gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
 {
   const size_t N = A->size1;
   size_t i, j;

   /* initial Cholesky */
   int stat_chol = gsl_linalg_cholesky_decomp(A);

   if(stat_chol == GSL_SUCCESS)
   {
     /* calculate D from diagonal part of initial Cholesky */
     for(i = 0; i < N; ++i)
     {
       const double C_ii = gsl_matrix_get(A, i, i);
       gsl_vector_set(D, i, C_ii*C_ii);
     }

     /* multiply initial Cholesky by 1/sqrt(D) on the right */
     for(i = 0; i < N; ++i)
     {
       for(j = 0; j < N; ++j)
       {
         gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
       }
     }

     /* Because the initial Cholesky contained both L and transpose(L),
        the result of the multiplication is not symmetric anymore;
        but the lower triangle _is_ correct. Therefore we reflect
        it to the upper triangle and declare victory.
        */
     for(i = 0; i < N; ++i)
       for(j = i + 1; j < N; ++j)
         gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
   }

   return stat_chol;
 }
	/* Cholesky Decomposition
	*
	* Copyright (C) 2000 Thomas Walter
	*
	* 03 May 2000: Modified for GSL by Brian Gough
	* 29 Jul 2005: Additions by Gerard Jungman
	* Copyright (C) 2000,2001, 2002, 2003, 2005 Brian Gough, Gerard Jungman
	*
	* This is free software; you can redistribute it and/or modify it
	* under the terms of the GNU General Public License as published by the
	* Free Software Foundation; either version 2, or (at your option) any
	* later version.
	*
	* This source 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 General Public License
	* for more details.
	*/

	/*
	* Cholesky decomposition of a symmetrix positive definite matrix.
	* This is useful to solve the matrix arising in
	* periodic cubic splines
	* approximating splines
	*
	* This algorithm does:
	* A = L * L'
	* with
	* L := lower left triangle matrix
	* L' := the transposed form of L.
	*
	*/

	#include <config.h>

	#include <gsl/gsl_math.h>
	#include <gsl/gsl_errno.h>
	#include <gsl/gsl_vector.h>
	#include <gsl/gsl_matrix.h>
	#include <gsl/gsl_blas.h>
	#include <gsl/gsl_linalg.h>

	static inline
	double
	quiet_sqrt (double x)
	/* avoids runtime error, for checking matrix for positive definiteness */
	{
	return (x >= 0) ? sqrt(x) : GSL_NAN;
	}

	int
	gsl_linalg_cholesky_decomp (gsl_matrix * A)
	{
	const size_t M = A->size1;
	const size_t N = A->size2;

	if (M != N)
	{
	GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
	}
	else
	{
	size_t i,j,k;
	int status = 0;

	/* Do the first 2 rows explicitly. It is simple, and faster. And
	* one can return if the matrix has only 1 or 2 rows.
	*/

	double A_00 = gsl_matrix_get (A, 0, 0);

	double L_00 = quiet_sqrt(A_00);

	if (A_00 <= 0)
	{
	status = GSL_EDOM ;
	}

	gsl_matrix_set (A, 0, 0, L_00);

	if (M > 1)
	{
	double A_10 = gsl_matrix_get (A, 1, 0);
	double A_11 = gsl_matrix_get (A, 1, 1);

	double L_10 = A_10 / L_00;
	double diag = A_11 - L_10 * L_10;
	double L_11 = quiet_sqrt(diag);

	if (diag <= 0)
	{
	status = GSL_EDOM;
	}

	gsl_matrix_set (A, 1, 0, L_10);
	gsl_matrix_set (A, 1, 1, L_11);
	}

	for (k = 2; k < M; k++)
	{
	double A_kk = gsl_matrix_get (A, k, k);

	for (i = 0; i < k; i++)
	{
	double sum = 0;

	double A_ki = gsl_matrix_get (A, k, i);
	double A_ii = gsl_matrix_get (A, i, i);

	gsl_vector_view ci = gsl_matrix_row (A, i);
	gsl_vector_view ck = gsl_matrix_row (A, k);

	if (i > 0) {
	gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
	gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);

	gsl_blas_ddot (&di.vector, &dk.vector, &sum);
	}

	A_ki = (A_ki - sum) / A_ii;
	gsl_matrix_set (A, k, i, A_ki);
	}

	{
	gsl_vector_view ck = gsl_matrix_row (A, k);
	gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);

	double sum = gsl_blas_dnrm2 (&dk.vector);
	double diag = A_kk - sum * sum;

	double L_kk = quiet_sqrt(diag);

	if (diag <= 0)
	{
	status = GSL_EDOM;
	}

	gsl_matrix_set (A, k, k, L_kk);
	}
	}

	/* Now copy the transposed lower triangle to the upper triangle,
	* the diagonal is common.
	*/

	for (i = 1; i < M; i++)
	{
	for (j = 0; j < i; j++)
	{
	double A_ij = gsl_matrix_get (A, i, j);
	gsl_matrix_set (A, j, i, A_ij);
	}
	}

	if (status == GSL_EDOM)
	{
	GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
	}

	return GSL_SUCCESS;
	}
	}


	int
	gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
	const gsl_vector * b,
	gsl_vector * x)
	{
	if (LLT->size1 != LLT->size2)
	{
	GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
	}
	else if (LLT->size1 != b->size)
	{
	GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
	}
	else if (LLT->size2 != x->size)
	{
	GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
	}
	else
	{
	/* Copy x <- b */

	gsl_vector_memcpy (x, b);

	/* Solve for c using forward-substitution, L c = b */

	gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

	/* Perform back-substitution, U x = c */

	gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);


	return GSL_SUCCESS;
	}
	}

	int
	gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
	gsl_vector * x)
	{
	if (LLT->size1 != LLT->size2)
	{
	GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
	}
	else if (LLT->size2 != x->size)
	{
	GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
	}
	else
	{
	/* Solve for c using forward-substitution, L c = b */

	gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

	/* Perform back-substitution, U x = c */

	gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

	return GSL_SUCCESS;
	}
	}


	int
	gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
	{
	const size_t N = A->size1;
	size_t i, j;

	/* initial Cholesky */
	int stat_chol = gsl_linalg_cholesky_decomp(A);

	if(stat_chol == GSL_SUCCESS)
	{
	/* calculate D from diagonal part of initial Cholesky */
	for(i = 0; i < N; ++i)
	{
	const double C_ii = gsl_matrix_get(A, i, i);
	gsl_vector_set(D, i, C_ii*C_ii);
	}

	/* multiply initial Cholesky by 1/sqrt(D) on the right */
	for(i = 0; i < N; ++i)
	{
	for(j = 0; j < N; ++j)
	{
	gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
	}
	}

	/* Because the initial Cholesky contained both L and transpose(L),
	the result of the multiplication is not symmetric anymore;
	but the lower triangle _is_ correct. Therefore we reflect
	it to the upper triangle and declare victory.
	*/
	for(i = 0; i < N; ++i)
	for(j = i + 1; j < N; ++j)
	gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
	}

	return stat_chol;
	}