home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
World_Of_Computer_Software-02-385-Vol-1of3.iso
/
x
/
xlmath10.zip
/
XLMATH.C
< prev
next >
Wrap
C/C++ Source or Header
|
1992-03-14
|
19KB
|
663 lines
/*********************************************************************
XLMATH.C
========
This module demonstrates the basic structure of a Dynamic Link Library
that can be used with Microsoft Excel.
Author: Roy Kari
Dept. of Chemistry
Laurentian University
Sudbury, Ont.
Canada P3E 2C6
(705) 675-1151
Internet: "ROY@NICKEL.LAURENTIAN.CA"
********************************************************************/
/* --------------------------< Include files >--------------------- */
#include <windows.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "xlmopti.h"
#include "xlmath.h"
#include "xlmutil.h"
#include "xlmcurve.h"
/* --------------------<forward references>------------------------ */
static WORD PASCAL BinBucket(LPFP lpValues, double dMin, double dMax);
static int PASCAL Jacobi(LPLPREAL hmat, LPLPREAL eivec,
LPREAL eival, int n);
static void PASCAL SortVectors(LPLPREAL eivec, LPREAL eival, int n);
/* ----------------<Frequency utility routine>---------------------- */
/*********************************************************************
BinBucket()
==============
This is an EXCEL math routine that calculates the number of occurences
greater than the lower bound and less than or equal to an upperbound.
This routine is used internally by Frequency() but can also be called
from EXCEL.
********************************************************************/
static WORD PASCAL BinBucket(LPFP lpValues, double dMin, double dMax)
{
WORD wValueRow;
WORD wBinCount = 0;
// no memory allocated so don't bother checking initialization
for (wValueRow = 0; wValueRow < lpValues->wRows; ++wValueRow)
{
if (lpValues->Data[wValueRow] > dMin &&
lpValues->Data[wValueRow] <= dMax)
{
++wBinCount;
}
}
return wBinCount;
}
/*********************************************************************
Frequency()
==============
An Excel DLL math routine to calculate the freqency of occurences.
Each value in the bin block represents all values from it down to
the previous value. The first bin value represents the number of values
less than or equal to itself. The last bin value represents the number
of values greater than itself. The number of bin values returned is
always one greater than the number of bin values.
********************************************************************/
LPFP PASCAL FAR _export Frequency(LPFP lpValues, LPFP lpBrackets)
{
WORD wValueRow, wBracketRow;
LPFP lpXLKInt;
// Init XL transfer array 1 row greater
if ( (lpXLKInt =InitRetBuff(lpBrackets->wRows+1, lpBrackets->wCols)) == NULL)
return (LPFP) NULL;
// zero the bin counts
for (wBracketRow = 0; wBracketRow <= lpBrackets->wRows; ++wBracketRow)
{
lpXLKInt->Data[wBracketRow] = 0;
}
// check first and last limits
for (wValueRow = 0; wValueRow < lpValues->wRows; ++wValueRow)
{
// value is <= first bin?
if ( lpValues->Data[wValueRow] <= lpBrackets->Data[0] )
{
++lpXLKInt->Data[0];
}
// value is > last bin?
else
if ( lpValues->Data[wValueRow] >
lpBrackets->Data[lpBrackets->wRows-1])
{
++lpXLKInt->Data[lpBrackets->wRows];
}
}
// check remaining bins if value > bin[-1] and <= bin[0]
for (wBracketRow = 1; wBracketRow < lpBrackets->wRows; ++wBracketRow)
{
lpXLKInt->Data[wBracketRow] = (double)BinBucket(lpValues,
lpBrackets->Data[wBracketRow-1],lpBrackets->Data[wBracketRow]);
}
return ((LPFP)lpXLKInt); // return as a long pointer
}
/* ----------------<Diagonalize utility routines>------------------ */
/********************************************************************
Jacobi()
========
This routine performs a diagonalization of a real symmetric matrix.
Only the the top half need be defined. The routine is adapted from
a FORTRAN routine of the same name. It is long and probably
should be broken into manageable portions but since it is
the biggest single time user in this application, I have opted
for speed over clarity.
*******************************************************************/
/* -------------------------< Defines >---------------------------- */
#define SQRTWO 0.707106781
#define ACCRCY DBL_EPSILON*1.e02 /* accuracy of diag */
#define MACHACC DBL_EPSILON /* machine precision */
static int PASCAL Jacobi(LPLPREAL hmat, LPLPREAL eivec,
LPREAL eival, int n)
{
int i, j; /* loop counters */
int converged; /* convergence control */
int jcol, irow; /* loop counters */
int it_cnt; /* iteration counter */
double htop, dstop, thrsh, d; /* stopping criteria etc */
double fd, avgf; /* temporary variable */
double hii, hjj, hij; /* elements of hmat */
double c, s, t, u, dsu, dcu; /* cos, sin, tan etc */
it_cnt = 0;
/* set eigenvector matrix to unit matrix */
for( i = 0; i < n; i++)
{
/* set eivec to unit matrix */
for (j = 0; j < n; j++)
{
eivec[i][j] = ( i==j ) ? 1.0: 0.0;
}
}
/* find max element of abs(hmat) */
htop = 0.0;
for(i = 0; i < n; i++)
{
for (j = i; j < n; j++)
{
htop = (htop > fabs(hmat[i][j])) ? htop : fabs(hmat[i][j]);
}
}
/* if only one element then done */
if (n == 1)
{
eival[0] = hmat[0][0];
return(0);
}
/* find the norm of the function and calculate stopping criteria */
avgf = (n-1)*n*0.55;
d = 0.0;
for (i = 1; i < n; i++)
{
for(j = i;j < n; j++)
{
s=hmat[i-1][j]/htop;
d = s*s + d;
}
}
dstop = d*ACCRCY;
thrsh = sqrt(d/avgf)*htop;
while (TRUE)
{
/* MAIN LOOP */
converged = TRUE;
/* loop over columns */
for (jcol=1; jcol<n; jcol++)
{
/* loop over rows */
for (irow=0; irow<jcol; irow++)
{
hij = hmat[irow][jcol];
if(fabs(hij) > thrsh )
{
/* executed for off-diagonal element > threshold */
hii = hmat[irow][irow];
hjj = hmat[jcol][jcol];
s = hjj - hii;
if (fabs(hij) > MACHACC*fabs(s))
{
/* executed for rotations greater then rounding error
and any element nulls convergence
*/
converged = FALSE;
if((double)(MACHACC*0.1)*fabs(hij) >= fabs(s))
{
/*
* If the rotation is close to 45 degrees,
* sin and cos to 1/root2.
* else calculate both sin & cos
*/
s = SQRTWO;
c = s;
}
else
{
t = hij/s;
s = 0.25/sqrt(0.25 + t*t);
c = sqrt( s+0.5);
s = 2.*t*s/c;
}
/*
* calculate new elements of hmat
*/
for (i=0; i<=irow; i++)
{
t = hmat[i][irow];
u = hmat[i][jcol];
dsu = s*u;
hmat[i][irow] = c*t-dsu;
dcu = c*u;
hmat[i][jcol] = s*t + dcu;
}
if( (irow+2) <= jcol )
{
for (i = irow+2; i <= jcol; i++)
{
t = hmat[i-1][jcol];
u = hmat[irow][i-1];
dsu = c*t;
hmat[i-1][jcol] = s*u + dsu;
dcu = s*t;
hmat[irow][i-1] = c*u - dcu;
}
}
dsu = c*hjj;
hmat[jcol][jcol] = s*hij + dsu;
dcu = s*hjj;
fd = c*hij - dcu;
dsu = s*fd;
hmat[irow][irow] = c* hmat[irow][irow] - dsu;
for (j=jcol; j<n; j++)
{
t = hmat[irow][j];
u = hmat[jcol][j];
dsu = s*u;
dcu = c*u;
hmat[irow][j] = c*t - dsu;
hmat[jcol][j] = s*t + dcu;
}
/*
* begin calculation of eigenvectors
*/
for (i=0; i<n; i++)
{
t = eivec[i][irow];
dsu = eivec[i][jcol]*s;
eivec[i][irow] = c*t - dsu;
dcu = eivec[i][jcol] * c;
eivec[i][jcol] = s*t + dcu;
}
/* recalculate the norm d and compare with dstop */
it_cnt += 1;
s = hij/htop;
d = d-s*s;
if( d < dstop )
{
d = 0.0;
for (i = 1;i < n; i++)
{
for(j = i;j < n; j++)
{
s=hmat[i-1][j]/htop;
d = s*s + d;
}
}
dstop = d*ACCRCY;
}
thrsh = sqrt(d/avgf)*htop;
} /* end rounding if */
} /* end threshold if */
} /* end row loop */
} /* end column loop */
if (converged)
break; /* end infinite while */
}
/* copy diagonal elements to eival */
for (i = 0; i < n; i++)
{
eival[i] = hmat[i][i];
}
/* get in beta order */
SortVectors(eivec, eival, n);
return(it_cnt);
}
/********************************************************************
*
* Name SortVectors -- sort a matrix
*
* Synopsis SortVectors(eivec, vec, n);
*
* double eivec[][] - eigenvector matrix
* (sorted to correpond to eivec)
* double vec[] - eigenvalue matrix (determines sort order)
* int n - no. of atoms
*
* Description (Shell) sort matrix eivec[][] in increasing order of
* values in vec[].
*
* Returns eivec[][] and eival[]
*
********************************************************************/
static void PASCAL SortVectors(LPLPREAL eivec, LPREAL eival, int n)
{
int gap;
register int i, j, k;
double temp;
for (gap = n/2; gap > 0; gap /=2) {
for ( i=gap; i<n; i++) {
for ( j=i-gap; j>=0 && eival[j] < eival[j+gap]; j-=gap) {
temp = eival[j]; /* exchange vec values */
eival[j] = eival[j+gap];
eival[j+gap] = temp;
for (k=0; k<n; k++) { /* exchange corresponding eivec */
temp = eivec[k][j]; /* columns */
eivec[k][j] = eivec[k][j+gap];
eivec[k][j+gap] = temp;
}
}
}
}
}
/********************************************************************
Diagonalize()
=============
Diagonalize a real symmetric matrix. Eigenvalues are returned
*******************************************************************/
LPFP PASCAL FAR _export Diagonalize(LPFP lpHmat)
{
LPFP lpXLKInt;
LPLPREAL lplpHmat; // pointers to EXCEL allocated data array
LPLPREAL lplpEivec; // handle to pointers
LPREAL lpEival; // pointer to eigenvalue array
WORD id; // index
int nIter; // no iterations in Jacobi
WORD wRows = lpHmat->wRows; // saves typing
BOOL Success = TRUE;
// error if not square array
if (wRows != lpHmat->wCols)
{
ErrorHandler(XLM_NOT_SQUARE);
return (LPFP)NULL;
}
// init XLK 1 row larger to store eigenvalues in last row
if ((lpXLKInt =InitRetBuff(wRows+1, wRows)) == NULL)
return (LPFP)NULL;
// create pointers for XLK internal array which is
// used to store vectors + eigenvalues in bottom row
if ((lplpEivec = InitPointers(&lpXLKInt->Data[0], wRows+1, wRows))
== NULL)
{
Success = FALSE;
goto Error3;
}
// create pointers for Hmat
if ((lplpHmat = InitPointers(&lpHmat->Data[0], wRows, wRows))
== NULL)
{
Success = FALSE;
goto Error2;
}
// allocate a Global 1-D array for eigenvalues
if ((lpEival = (LPREAL)GetMem( sizeof(double)*wRows )) == NULL)
{
Success = FALSE;
goto Error1;
}
// diagonalize the matrix
nIter = Jacobi(lplpHmat, lplpEivec, lpEival, wRows);
// copy the eigenvalues
for (id = 0; id < wRows; id++)
{
lplpEivec[wRows][id] = lpEival[id];
}
// free pointers and arrays
FreeMem((LPVOID)(LPVOID)lpEival);
Error1:
MemFreePtr((LPVOID)lplpEivec);
Error2:
MemFreePtr((LPVOID)lplpHmat);
Error3:
if (!Success)
{
return ((LPFP)NULL);
}
// return pointer to XL map to EXCEL
return ( (LPFP) lpXLKInt);
}
/********************************************************************
PolyCurveFit()
=============
This function is a generalized least squares curve fitting function.
It fits a polynomial with linear coefficients to a dependent -
independent variable set of data
*******************************************************************/
LPFP PASCAL FAR _export PolyCurveFit(LPFP lpIndVar, LPFP lpDepVar, unsigned int Order)
{
#define nXLKCols 3 // placed in a NumObs x nXLKCols array
LPFP lpXLKInt;
LPREAL lpYest, lpResid, lpPolycoef, lpCoefsig;
static double See, Rsqrval, Rval, dferror;
WORD wRows = lpIndVar->wRows;
WORD wCols = lpIndVar->wCols;
WORD NumObs = wRows > wCols ? wRows : wCols;
WORD wMaxRows = 2*(Order+1) + 4;
BOOL Success = TRUE;
WORD i, wIndex;
if (wMaxRows < wRows)
wMaxRows = wRows;
// init XLK & zero array data
if ((lpXLKInt = InitRetBuff(wMaxRows, nXLKCols)) == NULL)
{
return (LPFP)NULL;
}
for (i=0; i<wMaxRows*nXLKCols; i++)
lpXLKInt->Data[i] = 0.0;
// memory for estimated Y values
lpYest = (LPREAL)GetMem( NumObs*sizeof(double));
if (lpYest == NULL)
{
Success = FALSE;
goto Error4;
}
// memory for residuals
if ((lpResid = (LPREAL)GetMem( NumObs*sizeof(double))) == NULL)
{
Success = FALSE;
goto Error3;
}
//memory for coefficients of fitted polynomial
if ((lpPolycoef = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
{
Success = FALSE;
goto Error2;
}
// memory for standard errors of coefficient estimates
if ((lpCoefsig = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
{
Success = FALSE;
goto Error1;
}
Success = _PolyCurveFit((LPREAL)&lpIndVar->Data[0],
(LPREAL)&lpDepVar->Data[0], (int)NumObs, (int)Order, lpPolycoef,
lpYest, lpResid, (NPREAL)&See, lpCoefsig, &Rsqrval, &Rval, &dferror);
if (Success == TRUE)
{
// remaining variables stuck in XLKInt directly below polycoef
wIndex = 2*(Order+1);
lpXLKInt->Data[(wIndex)*nXLKCols + 2] = See;
lpXLKInt->Data[(wIndex+1)*nXLKCols + 2] = Rsqrval;
lpXLKInt->Data[(wIndex+2)*nXLKCols + 2] = Rval;
lpXLKInt->Data[(wIndex+3)*nXLKCols + 2] = dferror;
// copy data into XLKInt
for (i=0; i < NumObs; i++)
{
wIndex = i*nXLKCols;
lpXLKInt->Data[wIndex] = lpYest[i]; // Y estimates
lpXLKInt->Data[wIndex+1] = lpResid[i]; // residuals
if ( i <= Order )
{
lpXLKInt->Data[wIndex+2] = lpPolycoef[i]; //
wIndex = (i+Order+1)*nXLKCols;
lpXLKInt->Data[wIndex+2] = lpCoefsig[i]; //
}
}
}
FreeMem((LPVOID)lpCoefsig);
Error1:
FreeMem((LPVOID)lpPolycoef);
Error2:
FreeMem((LPVOID)lpResid);
Error3:
FreeMem((LPVOID)lpYest);
Error4:
if (!Success)
{
return ((LPFP) NULL);
}
return (lpXLKInt);
}
#undef nXLKCols
/********************************************************************
CubicSplines()
=============
This function will fit a set of cubic spline polynomial equations
to a discrete set of data.
*******************************************************************/
LPFP PASCAL FAR _export CubicSplines(LPFP lpIndVar, LPFP lpDepVar)
{
#define nXLKCols 4 // placed in a wRows x 4 array
LPFP lpXLKInt;
WORD wRows = lpIndVar->wRows;
WORD wCols = lpIndVar->wCols;
BOOL Success = TRUE;
WORD i, wIndex;
// init XLK & zero array data
if ((lpXLKInt = InitRetBuff(wRows-1, nXLKCols)) == NULL)
return (LPFP)NULL;
Success = _CubicSplines((LPREAL)&lpIndVar->Data[0],
(LPREAL)&lpDepVar->Data[0],
(int)wRows-1,
(LPREAL)&lpXLKInt->Data[0]);
if (!Success)
{
return ((LPFP)NULL);
}
return (lpXLKInt);
}
#undef nXLKCols
/********************************************************************
CalcSpline()
=============
This function calculate the cubic spline interpolation of an
y-value given an x-value and the cubic spline coefficients.
*******************************************************************/
LPREAL PASCAL FAR _export CalcSpline(LPFP lpIndVar, LPFP lpCoef, LPREAL lpX)
{
LPREAL lpXLEInt;
if ((lpXLEInt = (LPREAL)InitRetBuff(1,1)) == NULL)
return (NULL);
_CalcSpline((LPREAL)&lpIndVar->Data[0],
(LPREAL)&lpCoef->Data[0],
(int)(lpIndVar->wRows - 1),
(LPREAL)lpX,
(LPREAL)lpXLEInt);
return (lpXLEInt);
}
/********************************************************************
SmoothSG()
=============
This function uses the Savitsky - Golay algorithm to reduce the
noise in a sampled data set.
*******************************************************************/
LPFP PASCAL FAR _export SmoothSG(LPFP lpData, unsigned int wSmoothNum,
unsigned int wDerivNum)
{
LPFP lpXLKInt;
WORD wRows = lpData->wRows;
WORD wCols = lpData->wCols;
WORD wNumDat;
// init XLK & zero array data
if ((lpXLKInt = InitRetBuff(wRows, wCols)) == NULL)
return (LPFP)NULL;
// check limits
if (wSmoothNum < 1 || wSmoothNum > 5)
return (LPFP)NULL;
if (wDerivNum > 2)
return (LPFP)NULL;
// set NumDat
wNumDat = (wRows > wCols) ? wRows: wCols;
// call smoothing routine
_DataSmoothSg((LPREAL)&lpData->Data[0],
(int)wNumDat,
(int)wSmoothNum,
(int)wDerivNum,
(LPREAL)&lpXLKInt->Data[0]);
return (lpXLKInt);
}
/********************************************************************
SmoothWeights()
=============
This function smooths data via a weighted average
*******************************************************************/
LPFP PASCAL FAR _export SmoothWeights(LPFP lpData, LPFP lpWeights, LPREAL lpDivisor)
{
LPFP lpXLKInt;
WORD wRows = lpData->wRows;
WORD wCols = lpData->wCols;
WORD wNumDat;
WORD wSmoothNum;
// init XLK & zero array data
if ((lpXLKInt = InitRetBuff(wRows, wCols)) == NULL)
return (LPFP)NULL;
wNumDat = (wRows > wCols) ? wRows : wCols;
wSmoothNum = (lpWeights->wRows > lpWeights->wCols) ? lpWeights->wRows :
lpWeights->wCols;
_DataSmoothWeights((LPREAL)&lpData->Data[0],
(int)wNumDat,
(int)wSmoothNum,
(LPREAL)&lpWeights->Data[0],
(LPREAL) lpDivisor,
(LPREAL)&lpXLKInt->Data[0]);
return (lpXLKInt);
}