dcond.cc

来自「一个统计模块的C++源程序！应该能用到的！」· CC 代码 · 共 72 行

/*
Name          dcond - Improves the accuracy of dsweep when the diagonal 
                      elements of a are not of the same order of magnitude.

Usage         #include "usual.h"
              #include "tools.h"
              void dcond(REAL *a, INTEGER n, REAL *s,INTEGER isw)

Prototype in  tools.h

Description   a is a symmetric, positive definite, n by n matrix stored 
              columnwise with no unused space and first element in a[0]; i.e.,
              for (j=1; j<=n; j++) for (i=1; i<=n; i++) aij=a[n*(j-1)+(i-1)];
              will traverse the matrix with aij being the element in the 
              i-th row and j-th column.  s is a work vector of length n.
              The intended calling sequence is:
                dcond(a,n,s,0);
                ier=dsweep(a,n,eps);
                dcond(a,n,s,1);
              The usage
                dcond(a,n,s,0);
              will return a correlation matrix in a and the standard errors 
              in s when a is a variance-covariance matrix.

Remark        dscond.c is dcond.f as translated by f2c (version of 3 February 
              1990  3:36:42) with some reworking to remove dependence on f2c 
              libraries.

Return value  None.

Functions     Library: (none)
called        Sublib:  (none)
*/

#include "usual.h"
#include "tools.h"

void dcond(REAL *a, INTEGER n, REAL *s, INTEGER isw)
{
    /* System generated locals */
    INTEGER a_dim1, a_offset, i_1, i_2;

    /* Local variables */
    INTEGER i, j;

    /* Parameter adjustments */
    a_dim1 = n;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --s;

    /* Function Body */
    if (isw == 1) {
        goto L20;
    }
    i_1 = n;
    for (i = 1; i <= i_1; ++i) {
        s[i] = sqrt(a[i + i * a_dim1]);
    }
L20:
    i_1 = n;
    for (i = 1; i <= i_1; ++i) {
        i_2 = n;
        for (j = 1; j <= i_2; ++j) {
            if (s[i] * s[j] != 0.) {
                a[i + j * a_dim1] /= s[i] * s[j];
            }
        }
    }
    return;
}