#include <math.h>
#include <stdio.h>
#include "NR.h"

Macros
#define	TINY 1.0e-20

Functions
void	ludcmp (double *a, int n, int index, double *d)

void	lubksb (double *a, int n, int index, double *b)

Macro Definition Documentation

◆ TINY

#define TINY 1.0e-20

Definition at line 6 of file luc.c.

Function Documentation

◆ lubksb()

void lubksb	(	double **	a,
		int	n,
		int *	index,
		double *	b
	)

Definition at line 90 of file luc.c.

                                                      {
  int i, ii = 0, ip, j;
  double sum;
 
  for (i = 1; i <= n; i++) /* forward substitution. */
  {
    ip = index[i];
    sum = b[ip];
    b[ip] = b[i];
    if (ii) {
#ifdef _OPENMP
#pragma omp parallel for private(j) reduction(- : sum)
#endif
      for (j = ii; j <= i; j++) sum = sum - a[i][j] * b[j];
    } else if (sum)
      ii = i;
    b[i] = sum;
  }
 
  for (i = n; i >= 1; i--) /* back substitution. */
  {
    sum = b[i];
#ifdef _OPENMP
#pragma omp parallel for private(j) reduction(- : sum)
#endif
    for (j = i + 1; j <= n; j++) sum = sum - a[i][j] * b[j];
    b[i] = sum / a[i][i];
  }
}

Referenced by InvertMatrix().

◆ ludcmp()

void ludcmp	(	double **	a,
		int	n,
		int *	index,
		double *	d
	)

Definition at line 12 of file luc.c.

                                                      {
  int i, j, k;
  int imax = 1;
  double big, sum, dum;
  double *vv; /* vv stores the implicit scaling of each row */
 
  vv = dvector(1, n); /* No looping yet. */
  *d = 1.0;           /* loop over rows to get implicit scaling info. */
#ifdef _OPENMP
#pragma omp parallel for private(j, big)
#endif
  for (i = 1; i <= n; i++) {
    big = 0.0;
    for (j = 1; j <= n; j++)
      if (fabs(a[i][j]) > big) big = fabs(a[i][j]);
    if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
    vv[i] = 1.0 / big;
  }
 
  for (j = 1; j <= n; j++) /* outer loop of Crout's method */
  {
    for (i = 1; i < j; i++) {
      sum = a[i][j];
#ifdef _OPENMP
#pragma omp parallel for private(k) reduction(- : sum)
#endif
      for (k = 1; k < i; k++)  // OMPChk: parallelization may not help much
        sum = sum - a[i][k] * a[k][j];
      a[i][j] = sum;
    }
 
    big = 0.0; /* search for largest pivotal element. */
    for (i = j; i <= n; i++) {
      sum = a[i][j];
#ifdef _OPENMP
#pragma omp parallel for private(k) reduction(- : sum)
#endif
      for (k = 1; k < j; k++)  // OMPChk: parallelization may not help much
        sum = sum - a[i][k] * a[k][j];
      a[i][j] = sum;
 
      dum = vv[i] * fabs(sum);
      if (dum >= big) { /* is the figure of merit better */
        big = dum;      /* than the best? */
        imax = i;
      }
    }
 
    if (j != imax) {
#ifdef _OPENMP
#pragma omp parallel for private(k, dum)
#endif
      for (k = 1; k <= n; k++) {
        dum = a[imax][k];
        a[imax][k] = a[j][k];
        a[j][k] = dum;
      }
      *d = -(*d);       /* change the parity of d. */
      vv[imax] = vv[j]; /* interchange the scale factor. */
    }
    index[j] = imax;
    if (a[j][j] == 0.0) /* for some applications on singular */
      a[j][j] = TINY;   /* matrices, it is desirable to */
                        /* substitute TINY for zero. */
 
    if (j != n) /* finally, divide by the pivot element. */
    {
      dum = 1.0 / a[j][j];
#ifdef _OPENMP
#pragma omp parallel for private(i)  // OMPCheck: may not help much
#endif
      for (i = j + 1; i <= n; i++) a[i][j] = a[i][j] * dum;
    }
  }
 
  free_dvector(vv, 1, n);
}