The R Project SVN – R – Diff – /trunk/src/appl/uncmin.c

 /*
  *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1999-2017   The R Core Team
+ *  Copyright (C) 2003-2017   The R Foundation
  *  Copyright (C) 1997-1999   Saikat DebRoy
- *  Copyright (C) 1999-2015   The R Core Team
- *  Copyright (C) 2003-2010   The R Foundation
+ *
  *  This program 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 of the License, or
  *  (at your option) any later version.
 } /* lltslv */
 static void
 choldc(int nr, int n, double *a, double diagmx, double tol, double *addmax)
+{
-/* Find the perturbed l(l-transpose) [written ll+] decomposition
+/* Find the perturbed L(L-transpose) [written LL+] decomposition
  * of a+d, where d is a non-negative diagonal matrix added to a if
  * necessary to allow the cholesky decomposition to continue.
  * PARAMETERS :
  *	nr	     --> row dimension of matrix
  *	n	     --> dimension of problem
  *	a(n,n)	    <--> on entry: matrix for which to find perturbed
  *			      cholesky decomposition
- *			 on exit:  contains l of ll+ decomposition
+ *			 on exit:  contains L of LL+ decomposition
  *			 in lower triangular part and diagonal of "a"
  *	diagmx	     --> maximum diagonal element of "a"
  *	tol	     --> tolerance
  *	addmax	    <--	 maximum amount implicitly added to diagonal of "a"
  *			 in forming the cholesky decomposition of a+d
  *	internal variables
- *	aminl	 smallest element allowed on diagonal of l
+ *	aminl	 smallest element allowed on diagonal of L
  *	amnlsq	 =aminl**2
  *	offmax	 maximum off-diagonal element in column of a
  *	description
  *	the normal cholesky decomposition is performed.	 however, if at any
- *	point the algorithm would attempt to set l(i,i)=sqrt(temp)
+ *	point the algorithm would attempt to set L(i,i)=sqrt(temp)
- *	with temp < tol*diagmx, then l(i,i) is set to sqrt(tol*diagmx)
+ *	with temp < tol*diagmx, then L(i,i) is set to sqrt(tol*diagmx)
  *	instead.  this is equivalent to adding tol*diagmx-temp to a(i,i)
  */
     double tmp1, tmp2;
     int i, j, k;
     *addmax = 0.0;
     aminl = sqrt(diagmx * tol);
     amnlsq = aminl * aminl;
-    /*	form row i of l */
+    /*	form row i of L */
     for (i = 0; i < n; ++i) {
+	// A[i,j] := * || find i,j element of lower triangular matrix L
+	for (j = 0; j < i; ++j) {
+	    sum = 0.;
+	    for (k = 0; k < j; ++k)
+		sum += a[i + k * nr] * a[j + k * nr];
+	    a[i + j * nr] = (a[i + j * nr] - sum) / a[j + j * nr];
+	}
-	/*	find diagonal elements of l */
+	// A[i,i] := * || find diagonal elements of L
 	sum = 0.;
 	for (k = 0; k < i; ++k)
 	    sum += a[i + k * nr] * a[i + k * nr];
 	tmp1 = a[i + i * nr] - sum;
-	if (tmp1 >= amnlsq) {
+	if (tmp1 >= amnlsq) { // normal Cholesky
 	    a[i + i * nr] = sqrt(tmp1);
+	}
-	else {
+	else { // augment diagonal of L
 	    /*	find maximum off-diagonal element in row */
 	    offmax = 0.;
 	    for (j = 0; j < i; ++j) {
 		if(offmax < (tmp2 = fabs(a[i + j * nr])))
 		    offmax = tmp2;
 	    /* add to diagonal element to
 	     * allow cholesky decomposition to continue */
 	    a[i + i * nr] = sqrt(offmax);
 	    if(*addmax < (tmp2 = offmax - tmp1)) *addmax = tmp2;
+	}
-	/*	find i,j element of lower triangular matrix */
-	for (j = 0; j < i; ++j) {
-	    sum = 0.;
-	    for (k = 0; k < j; ++k)
-		sum += a[i + k * nr] * a[j + k * nr];
-	    a[i + j * nr] = (a[i + j * nr] - sum) / a[j + j * nr];
-	}
+    }
 } /* choldc */
 static void qraux1(int nr, int n, double *r, int i)
+{
 } /* secfac */
 static void
 chlhsn(int nr, int n, double *a, double epsm, double *sx, double *udiag)
+{
-/*	find the l(l-transpose) [written ll+] decomposition of the perturbed
+/*	find the l(l-transpose) [written LL+] decomposition of the perturbed
  *	model hessian matrix a+mu*i(where mu\0 and i is the identity matrix)
  *	which is safely positive definite.  if a is safely positive definite
  *	upon entry, then mu=0.
  * PARAMETERS :
  *	nr	     --> row dimension of matrix
  *	n	     --> dimension of problem
  *	a(n,n)	    <--> on entry; "a" is model hessian (only lower
  *			 triangular part and diagonal stored)
- *			 on exit:  a contains l of ll+ decomposition of
+ *			 on exit:  a contains l of LL+ decomposition of
  *			 perturbed model hessian in lower triangular
  *			 part and diagonal and contains hessian in upper
  *			 triangular part and udiag
  *	epsm	     --> machine epsilon
  *	sx(n)	     --> diagonal scaling matrix for x
  *	such that the diagonal of a:=a+mu*i is all positive
  *	with the ratio of its smallest to largest element on the
  *	order of sqrt(epsm).
  *	2. "a" undergoes a perturbed cholesky decomposition which
- *	results in an ll+ decomposition of a+d, where d is a
+ *	results in an LL+ decomposition of a+d, where d is a
  *	non-negative diagonal matrix which is implicitly added to
  *	"a" during the decomposition if "a" is not positive definite.
  *	"a" is retained and not changed during this process by
  *	copying l into the upper triangular part of "a" and the
  *	diagonal into udiag.  then the cholesky decomposition routine
     /*	step3
      *	if addmax=0, "a" was positive definite going into step 2,
-     *	the ll+ decomposition has been done, and we return.
+     *	the LL+ decomposition has been done, and we return.
      *	otherwise, addmax>0.  perturb "a" so that it is safely
-     *	diagonally dominant and find ll+ decomposition */
+     *	diagonally dominant and find LL+ decomposition */
     if (addmax > 0.0) {
 	/*	restore original "a" (lower triangular part and diagonal) */

The R Project SVN R

(root)/trunk/src/appl/uncmin.c – Rev 68947 → 72555