Witam serdecznie
Musze napisać program wyznaczający wartości własne macierzy. Znalazłem przykład,który podobno oblicza wartości i wektory własne, ale nie jestem w stanie go odpalić, ponieważ pokazuje mi pewne błędy, których niestety z moją "wiedzą" (raczej niewiedzą) z c++ nie jestem w stanie naprawić. Program składa sie z 4 plików (1 cpp główny, drugi cpp z definicją funkcji jacobiego która oblicza te wartości, oraz z dwóch plików o rozszerzeniu.h ,które nie wiem co robią, bo sie nie znam na tym :/ )
Problem polega na tym ze podczas kompilacji dostaję ostrzeżenie:
[C++ Warning] ujacobi.cpp(20): W8058 Cannot create pre-compiled header: code in header
UJACOBI.cpp to plik zawierajacy definicje funkcji jacobiego
a przy próbie odpalenia programu wyskakuje mi:
[Linker Error] Unresolved external 'vminit()' referenced from D:\DOKUMENTY\STUDIA DOC\PSM\ZALICZENIE\TUJACOBI.OBJ
[Linker Error] Unresolved external 'vmalloc(void *, int, unsigned int, unsigned int)' referenced from D:\DOKUMENTY\STUDIA DOC\PSM\ZALICZENIE\TUJACOBI.OBJ
[Linker Error] Unresolved external 'vmfree(void *)' referenced from D:\DOKUMENTY\STUDIA DOC\PSM\ZALICZENIE\TUJACOBI.OBJ
TUJACOBI.cpp to plik głowny
pracuje w Builderze 6
Poniżej podam kody tych wszystkich plików.
Byłbym wdzięczny gdyby ktoś pomógł mi rozszyfrować te błędy, ponieważ jest to moja pierwsza styczność z programem składającym sie z wielu plików, więc pewnie popełniam jakiś głupi błąd. :|
Niestety na zajęciach robiliśmy proste programiki liczące pierwiastki funkcji kwadratowej, a na zaliczenie dostalem takie cudo. [glowa]
TUJACOBI.OBJ
/*********************************************************************
* Eigenvalues and eigenvectors of a real symmetric *
* square matrix by JACOBI'S METHOD *
* ------------------------------------------------------------------ *
* Uses function Jacobi of module ujacobi.cpp. *
* *
* C++ version by J-P Moreau *
* ------------------------------------------------------------------ *
* SAMPLE RUN: *
* Input file (elpro4.dat): *
* *
* 4 *
* 4 -2 -1 0 *
* 4 0 -1 *
* 4 -2 *
* 4 *
* *
* Output file (elpro.lst): *
* *
* --------------------------------------------------- *
* Eigenvalues and eigenvectors of a real, symmetric *
* square matrix by Jacobi's iterative method *
* (with dynamic allocations) *
* --------------------------------------------------- *
* *
* Dimension of matrix: 4 *
* *
* Symmetric matrix is: *
* *
* 4.000000 -2.000000 -1.000000 0.000000 *
* -2.000000 4.000000 0.000000 -1.000000 *
* -1.000000 0.000000 4.000000 -2.000000 *
* 0.000000 -1.000000 -2.000000 4.000000 *
* *
* Eigenvalues: *
* 7.00000000 *
* 1.00000000 *
* 3.00000000 *
* 5.00000000 *
* *
* Eigenvectors (in lines): *
* 1.000000 -1.000000 -1.000000 1.000000 *
* 1.000000 1.000000 1.000000 1.000000 *
* 1.000000 1.000000 -1.000000 -1.000000 *
* -1.000000 1.000000 -1.000000 1.000000 *
* *
* Number of rotations: 25 *
* --------------------------------------------------- *
* *
**********************************************************************
Version with dynamic allocations (to link with files basis_r.cpp and
vmblock.cpp). */
#include "ujacobi.cpp" //for function Jacobi()
void main() {
int dimen; // Dimension of square matrix
int nrot; // Number of matrix rotations
int i,j; // loop variables
double vmax; // normalization value
REAL **Mat; // given square symmetric matrix
REAL *Eigenvalues; // eigenvalues solution vector
REAL **Eigenvectors; // eigenvectors solution matrix
FILE *fp1, *fp2;
void *vmblock = NULL;
// open input/output files
fp1 =fopen("elpro4.dat","r"); // read mode
fp2 =fopen("elpro.lst","w"); // write mode
//read size of matrix from input file
fscanf(fp1,"%d", &dimen);
//allocate vector and matrices
vmblock = vminit();
Mat = (REAL **) vmalloc(vmblock, MATRIX, dimen+1, dimen+1);
Eigenvalues = (REAL *) vmalloc(vmblock, VEKTOR, dimen+1, 0);
Eigenvectors = (REAL **) vmalloc(vmblock, MATRIX, dimen+1, dimen+1);
//read data from input file
for (i=1; i<=dimen; i++) // read only upper half triangle
for (j=i; j<=dimen; j++) {
fscanf(fp1,"%lf", &Mat[i][j]);
Mat[j][i]=Mat[i][j];
}
fclose(fp1);
//print title and input matrix
fprintf(fp2,"---------------------------------------------------\n");
fprintf(fp2," Eigenvalues and eigenvectors of a real, symmetric\n");
fprintf(fp2," square matrix by Jacobi's iterative method\n");
fprintf(fp2," (with dynamic allocations)\n");
fprintf(fp2,"---------------------------------------------------\n\n");
fprintf(fp2," Dimension of matrix: %d\n\n", dimen);
fprintf(fp2," Symmetric matrix is:\n");
for (i=1; i<=dimen; i++) {
for (j=1; j<=dimen; j++)
fprintf(fp2,"%10.6f", Mat[i][j]);
fprintf(fp2,"\n");
}
Jacobi(Mat,dimen,Eigenvalues,Eigenvectors,&nrot);
//print results to output file
fprintf(fp2,"\n\n Eigenvalues:\n");
for (i=1; i<=dimen; i++)
fprintf(fp2,"%14.8f\n", Eigenvalues[i]);
fprintf(fp2,"\n\n Eigenvectors (in lines):\n");
for (j=1; j<=dimen; j++) {
vmax=Eigenvectors[1][j];
for (i=2; i<=dimen; i++)
if(fabs(Eigenvectors[i][j])>fabs(vmax)) vmax=Eigenvectors[i][j];
for (i=1; i<=dimen; i++)
fprintf(fp2,"%10.6f", Eigenvectors[i][j]/vmax);
fprintf(fp2,"\n");
}
fprintf(fp2,"\n\n Number of rotations: %d\n", nrot);
fprintf(fp2,"---------------------------------------------------\n\n");
fclose(fp2);
printf("\n Program done.\n");
printf(" Results in file elpro.lst.\n\n");
vmfree(vmblock);
}
// end of file tujacobi.cpp
UJACOBI.cpp
#include "basis.h"
#include "vmblock.h"
/************************************************************
* This subroutine computes all eigenvalues and eigenvectors *
* of a real symmetric square matrix A(N,N). On output, ele- *
* ments of A above the diagonal are destroyed. D(N) returns *
* the eigenvalues of matrix A. V(N,N) contains, on output, *
* the eigenvectors of A by columns. THe normalization to *
* unity is made by main program before printing results. *
* NROT returns the number of Jacobi matrix rotations which *
* were required. *
* --------------------------------------------------------- *
* Ref.:"NUMERICAL RECIPES IN FORTRAN, Cambridge University *
* Press, 1986, chap. 11, pages 346-348" [BIBLI 08]. *
* *
* C++ version by J-P Moreau, Paris. *
************************************************************/
void Jacobi(REAL **A,int N,REAL *D, REAL **V, int *NROT) {
REAL *B, *Z;
double c,g,h,s,sm,t,tau,theta,tresh;
int i,j,ip,iq;
void *vmblock1 = NULL;
//allocate vectors B, Z
vmblock1 = vminit();
B = (REAL *) vmalloc(vmblock1, VEKTOR, 100, 0);
Z = (REAL *) vmalloc(vmblock1, VEKTOR, 100, 0);
for (ip=1; ip<=N; ip++) { //initialize V to identity matrix
for (iq=1; iq<=N; iq++) V[ip][iq]=0;
V[ip][ip]=1;
}
for (ip=1; ip<=N; ip++) {
B[ip]=A[ip][ip];
D[ip]=B[ip];
Z[ip]=0;
}
*NROT=0;
for (i=1; i<=50; i++) {
sm=0;
for (ip=1; ip<N; ip++) //sum off-diagonal elements
for (iq=ip+1; iq<=N; iq++)
sm=sm+fabs(A[ip][iq]);
if (sm==0) {
vmfree(vmblock1);
return; //normal return
}
if (i < 4)
tresh=0.2*sm*sm;
else
tresh=0;
for (ip=1; ip<N; ip++) {
for (iq=ip+1; iq<=N; iq++) {
g=100*fabs(A[ip][iq]);
// after 4 sweeps, skip the rotation if the off-diagonal element is small
if ((i > 4) && (fabs(D[ip])+g == fabs(D[ip])) && (fabs(D[iq])+g == fabs(D[iq])))
A[ip][iq]=0;
else if (fabs(A[ip][iq]) > tresh) {
h=D[iq]-D[ip];
if (fabs(h)+g == fabs(h))
t=A[ip][iq]/h;
else {
theta=0.5*h/A[ip][iq];
t=1/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0) t=-t;
}
c=1.0/sqrt(1.0+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*A[ip][iq];
Z[ip] -= h;
Z[iq] += h;
D[ip] -= h;
D[iq] += h;
A[ip][iq]=0;
for (j=1; j<ip; j++) {
g=A[j][ip];
h=A[j][iq];
A[j][ip] = g-s*(h+g*tau);
A[j][iq] = h+s*(g-h*tau);
}
for (j=ip+1; j<iq; j++) {
g=A[ip][j];
h=A[j][iq];
A[ip][j] = g-s*(h+g*tau);
A[j][iq] = h+s*(g-h*tau);
}
for (j=iq+1; j<=N; j++) {
g=A[ip][j];
h=A[iq][j];
A[ip][j] = g-s*(h+g*tau);
A[iq][j] = h+s*(g-h*tau);
}
for (j=1; j<=N; j++) {
g=V[j][ip];
h=V[j][iq];
V[j][ip] = g-s*(h+g*tau);
V[j][iq] = h+s*(g-h*tau);
}
*NROT=*NROT+1;
} //end ((i.gt.4)...else if
} // main iq loop
} // main ip loop
for (ip=1; ip<=N; ip++) {
B[ip] += Z[ip];
D[ip]=B[ip];
Z[ip]=0;
}
} //end of main i loop
printf("\n 50 iterations !\n");
vmfree(vmblock1);
return; //too many iterations
}
// end of file ujacobi.cpp
basis.h
// Header file of basis_r.cpp
/***********************************************************************
* *
* Basic functions: Declarations file (with types and macros) *
* ---------------------------------------------------------- *
* *
* Programming language: ANSI C *
* Compiler: Turbo C 2.0 *
* Computer: IBM PS/2 70 with 80387 *
* Author: Juergen Dietel, Computer Center, RWTH Aachen *
* [BIBLI 11]. *
* Date: 9.30.1992 *
* *
***********************************************************************/
/***********************************************************************
* Make preparations in case this declarations file is included more *
* than once in a source code *
***********************************************************************/
#ifndef BASIS_H_INCLUDED
#define BASIS_H_INCLUDED
/***********************************************************************
* Include a few other declarations file here. *
* We include the standard declarations files here mainly so that only *
* this file needs to be included in the C modules of the Numerical *
* Library in order to have access to all standard names and variables. *
* Moreover modifications for nonstandard compilers may then be limited *
* to this file only. *
***********************************************************************/
#ifdef __EMX__ /* emx 0.9b with GNU CC 2.7.2? */
#if __GNUC__ == 2 /* To prevent a compiler crash in */
#if __GNUC_MINOR__ == 7 /* `15\gax.c' while compiling with LDOUBLE, */
#define EMX09B /* by defining the macro _WITH_UNDERSCORE */
#define _WITH_UNDERSCORE /* the special variants of the mathematical */
#define fabsl _fabsl /* functions for `long double' in `math.h' */
#define sqrtl _sqrtl /* are activated and their use further */
#define powl _powl /* below in shape of the macros FABS, SQRT, */
#define sinl _sinl /* ... is prepared here; which seems quite */
#define cosl _cosl /* useful even without the affair with the */
#define expl _expl /* crashing compiler. */
#define logl _logl
#define atanl _atanl
#define acosl _acosl
#define coshl _coshl
#endif
#endif
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h> /* for fabs, sqrt, pow, exp, sin, cos, log, */
/* atan, acos */
#include <float.h> /* for DBL_EPSILON, DBL_MAX */
#include <string.h>
#ifdef sun /* for GNU CC under SunOS */
#include <unistd.h> /* for SEEK_END */
#else
#ifdef amigados /* for GNU CC on an Amiga */
#include <unistd.h> /* for SEEK_END */
#endif
#endif
#ifndef SEEK_END /* SEEK_END not yet defined (such as for */
/* GNU CC 2.1 for an i386 MS-DOS computer)? */
#endif
/***********************************************************************
* Predefine the desired precision for floating point arithmetic: *
* If the macro FLOAT has been defined, we compute in single precision *
* (data type float), for LDOUBLE we compute in maximal precision *
* (data type long double), otherwise we use double precision (data type*
* double). LDOUBLE only adds precision on those compilers for which *
* long double is different from double such as on Turbo C, but not on *
* QuickC compilers for example. *
* For input and output of floating numbers the macros LZS (length *
* prefix for scanf()) and LZP (length prefix for printf()) should be *
* used (obviously two different ones), since according to the ANSI C *
* Standard the output of double values should be done without the *
* length declaration "l" in contrast to the input where this "l" is *
* needed. *
* NOTE: If the compiler used does not know the macros FLT_MAX, *
* ==== LDBL_MAX, DBL_MAX or FLT_MAX_EXP, LDBL_MAX_EXP, DBL_MAX_EXP *
* (which can usually be found in float.h), the user must *
* insert suitable values at the locations marked with !!!!. *
***********************************************************************/
#ifdef FLOAT /* single precision ? ...............*/
typedef float REAL; /* Standard floating point type float*/
/*.IX{REAL}*/
typedef double LONG_REAL; /* more precise floating point type */
/*.IX{LONG\unt REAL}*/
#ifdef FLT_EPSILON /* ANSI C compiler ? ................*/
#define MACH_EPS (REAL)FLT_EPSILON
/*.IX{MACH\unt EPS}*/
#else /* not an ANSI C compiler ? .........*/
#define MACH_EPS mach_eps() /* machine constant .................*/
#endif
#ifdef FLT_MAX_EXP /* ANSI C compiler? .................*/
#define MAX_EXP FLT_MAX_EXP /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else /* not an ANSI C compiler ? .........*/
#define MAX_EXP 128 /* make adjustments here !!!! .......*/
#endif
#ifdef FLT_MAX /* ANSI C compiler? ................ */
#define POSMAX (REAL)FLT_MAX /* largest floating point number ... */
/*.IX{POS\unt MAX}*/
#else /* not an ANSI C compiler ? ........ */
#define POSMAX 1e38f /* adjust here !!!! ................ */
#endif
#ifdef FLT_MIN /* ANSI C compiler? .................*/
#define POSMIN (REAL)FLT_MIN /* smallest positive floating point */
/*.IX{POS\unt MIN}*/
/* number ...........................*/
#else /* not an ANSI C compiler ? .........*/
#define POSMIN posmin()
#endif
#define LZS "" /* length prefix for formatted */
/*.IX{LZS}*/
/* input of floating point numbers ..*/
#define LZP "" /* length prefix for formatted */
/*.IX{LZP}*/
/* output of floating point numbers .*/
#else
#ifdef LDOUBLE /* maximal precision ? ..............*/
typedef long double REAL; /* Standard floating point type */
/* long double */
/*.IX{REAL}*/
typedef long double LONG_REAL; /* "more precise" floating point type*/
/*.IX{LONG\unt REAL}*/
#define LONG_DOUBLE_USED
#ifdef LDBL_EPSILON /* ANSI C compiler? .................*/
#define MACH_EPS (REAL)LDBL_EPSILON
/*.IX{MACH\unt EPS}*/
#else /* not an ANSI C compiler ? .........*/
#define MACH_EPS mach_eps() /* machine constant .................*/
#endif
#ifdef LDBL_MAX_EXP /* ANSI C compiler? .................*/
#define MAX_EXP LDBL_MAX_EXP /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else /* not an ANSI C compiler ? .........*/
#define MAX_EXP 1023 /* adjust here !!!! .................*/
#endif
#ifdef LDBL_MAX /* ANSI C compiler? .................*/
#define POSMAX (REAL)LDBL_MAX /* largest floating point number ....*/
/*.IX{POS\unt MAX}*/
#else /* not an ANSI C compiler ? .........*/
#define POSMAX 1e100l /* adjust here !!!! .................*/
#endif
#ifdef LDBL_MIN /* ANSI C compiler? .................*/
#define POSMIN (REAL)LDBL_MIN /* smallest positive floating point */
/*.IX{POS\unt MIN}*/ /* number ...........................*/
#else /* not an ANSI C compiler ? ........*/
#define POSMIN posmin()
#endif
#define LZS "L" /* length prefix for formatted */
/*.IX{LZS}*/
/* input of floating point numbers ..*/
#define LZP "L" /* length prefix for formatted */
/*.IX{LZP}*/
/* output of floating point numbers .*/
#else /* double precision ? ...............*/
typedef double REAL; /* Standard floating point type */
/* double */
/*.IX{REAL}*/
typedef long double LONG_REAL; /* more precise floating point type */
/*.IX{LONG\unt REAL}*/
#ifdef DBL_EPSILON /* ANSI C compiler? .................*/
#define MACH_EPS (REAL)DBL_EPSILON
/*.IX{MACH\unt EPS}*/
#else /* not an ANSI C compiler ? .........*/
#define MACH_EPS mach_eps() /* machine constant .................*/
#endif
#ifdef DBL_MAX_EXP /* ANSI C compiler? .................*/
#define MAX_EXP DBL_MAX_EXP /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else /* not an ANSI C compiler ? .........*/
#define MAX_EXP 1023 /* adjust here !!!! .................*/
#endif
#ifdef DBL_MAX /* ANSI C compiler? .................*/
#define POSMAX (REAL)DBL_MAX /* largest floating point number ....*/
/*.IX{POS\unt MAX}*/
#else /* not an ANSI C compiler ? .........*/
#define POSMAX 1e100 /* adjust here !!!! .................*/
#endif
#ifdef DBL_MIN /* ANSI C compiler? .................*/
#ifdef __BORLANDC__
#if __BORLANDC__ <= 0x0200 /* Borland C++ 2.0 for DOS? .........*/
/* because the value */
#define POSMIN 2.2250738585072017E-308 /* 2.2250738585072014E-308 */
#else /* from `float.h' is */
/* regarded as zero! */
#define POSMIN DBL_MIN /* smallest positive floating point */
#endif /* number ...........................*/
#else
#define POSMIN DBL_MIN /* smallest positive floating point */
#endif
/*.IX{POS\unt MIN}*/ /* number ...........................*/
#else /* not an ANSI C compiler ? .........*/
#define POSMIN posmin()
#endif
#define LZS "l" /* length prefix for formatted */
/*.IX{LZS}*/
/* input of floating point numbers ..*/
#define LZP "" /* length prefix for formatted */
/*.IX{LZP}*/
/* output of floating point numbers .*/
#endif
#endif
/***********************************************************************
* declare several important data types *
***********************************************************************/
#undef boolean
#undef FALSE
#undef TRUE
typedef enum {FALSE, TRUE} boolean;
/*.IX{FALSE}*/
/*.IX{TRUE}*/
/*.IX{boolean}*/
/* function pointer types for approximation in chapter 8 .............*/
typedef REAL (*ansatzfnk) (int i, REAL x);
/*.IX{ansatzfnk}*/
typedef REAL (*approxfnk) (REAL c[], REAL x);
/*.IX{approxfnk}*/
typedef void (*ableitfnk) (REAL x, REAL c[], REAL *d);
/*.IX{ableitfnk}*/
/* Type of function, which evaluates the right hand side of an .......*/
/* explicit ordinary differential equation y' = f(x,y) .......*/
typedef REAL (*dglfnk)(REAL x, REAL y);
/*.IX{dglfnk}*/
/* Type of function, which evaluates the right hand side .............*/
/* of a system of differential equations y' = f(x,y) .............*/
typedef void (*dglsysfnk)(REAL x, REAL y[], REAL f[]);
/*.IX{dglsysfnk}*/
/* Type of function, which computes the boundary value r(ya, yb) .....*/
/* of a two-point first order boundary value problem .....*/
typedef void (*rndbedfnk)(REAL ya[], REAL yb[], REAL r[]);
/*.IX{rndbedfnk}*/
/* enumeration type for classification of error codes that are ......*/
/* returned from most of the functions that realize a numerical ......*/
/* method ......*/
typedef enum { KEIN_FEHLER, WARNUNG, UNBEKANNT, FATAL } fehler_t;
/*.IX{KEIN\unt FEHLER}*/
/*.IX{WARNUNG}*/
/*.IX{UNBEKANNT}*/
/*.IX{FATAL}*/
/*.IX{fehler\unt t}*/
typedef REAL abl_mat1[4][2]; /* used to evaluate splinefunctions */
/*.IX{abl\unt mat1}*/
typedef REAL abl_mat2[6][2]; /* in spliwert ......................*/
/*.IX{abl\unt mat2}*/
typedef REAL mat4x4[4][4]; /* type for bicubic splines */
/*--------------------------------------------------------------------*
* Type declarations by Albert Becker *
*--------------------------------------------------------------------*/
/* Real functions ...................................................*/
typedef REAL (* REALFCT) (REAL);
/*.IX{REALFCT}*/
/* Real multi-dimensional functions .................................*/
typedef int (* FNFCT) (int, REAL [], REAL []);
/*.IX{FNFCT}*/
/* Functions for finding the Jacobi matrix ..........................*/
typedef int (* JACOFCT) (int, REAL [], REAL * []);
/*.IX{JACOFCT}*/
/***********************************************************************
* define sevaral important macros *
* NOTE: Borland C++ offers floating point standard functions *
* suitable for long double type from version 3.0 on such as *
* expl() instead of exp(), sinl() instead of sin() etc. But *
* Borland C++ 3.0 does not seem to define a macro that lets *
* the user differentiate it from Borland C++ 2.0 and below, *
* the user is forced to do so himself: If using Borland C++ *
* 3.0 and long double the user should define the macro BC3 *
* before compiling. Only in this case will the new maximally *
* precise floating point functions be used. *
***********************************************************************/
#define BASIS 2 /* Basis of number representation */
/*.IX{BASIS}*/
#define EPSROOT epsroot() /* square root of MACH_EPS */
/*.IX{EPSROOT}*/
#define EPSQUAD epsquad() /* square of MACH_EPS */
/*.IX{EPSQUAD}*/
#define MAXROOT maxroot() /* square root of the largest */
/*.IX{MAXROOT}*/
/* floating point number */
#ifndef PI
#define PI pi() /* pi = 3.14... */
/*.IX{PI}*/
#endif
#define EXP_1 exp_1() /* e = 2.71... */
/*.IX{EXP\unt 1}*/
#define ZERO (REAL)0.0 /* declare names for recurring */
/*.IX{ZERO}*/
#define ONE (REAL)1.0 /* floating point numbers */
/*.IX{ONE}*/
#define TWO (REAL)2.0
/*.IX{TWO}*/
#define THREE (REAL)3.0
/*.IX{THREE}*/
#define FOUR (REAL)4.0
/*.IX{FOUR}*/
#define FIVE (REAL)5.0
/*.IX{FIVE}*/
#define SIX (REAL)6.0
/*.IX{SIX}*/
#define EIGHT (REAL)8.0
/*.IX{EIGHT}*/
#define NINE (REAL)9.0
/*.IX{NINE}*/
#define TEN (REAL)10.0
/*.IX{TEN}*/
#define HALF (REAL)0.5
/*.IX{HALF}*/
#ifdef __BORLANDC__
#if __BORLANDC__ >= 0x0400 /* a BC distinguishable from BC2 */
/* and with long double functions */
/* (at least Borland C++ 3.1 or */
/* Borland C++ 1.00 for OS/2)? */
#define BC3 /* define `BC3' automatically */
#endif
#endif
#ifdef _MSC_VER
#if _MSC_VER >= 0x0258 /* a MC distinguishable from QC2 */
/* and with long double functions */
/* (at least Microsoft C 6.00A)? */
#define MC6 /* define MC6 automatically */
#endif
#endif
#if defined(LDOUBLE) && /* Borland C++ 3.0 or */\
(defined(BC3) || defined(MC6) || /* Microsoft C 6.0 or */\
defined(EMX09B)) /* emx 0.9b (GCC272) with */
/* maximal precision? */
#define FABS(x) fabsl((x)) /* use the long double */
/*.IX{FABS}*/
#define SQRT(x) sqrtl((x)) /* versions of the basic */
/*.IX{SQRT}*/
#define POW(x, y) powl((x), (y)) /* floating point */
/*.IX{POW}*/
#define SIN(x) sinl((x)) /* functions */
/*.IX{SIN}*/
#define COS(x) cosl((x))
/*.IX{COS}*/
#define EXP(x) expl((x))
/*.IX{EXP}*/
#define LOG(x) logl((x))
/*.IX{LOG}*/
#define ATAN(x) atanl((x))
/*.IX{ATAN}*/
#define ACOS(x) acosl((x))
/*.IX{ACOS}*/
#define COSH(x) coshl((x))
/*.IX{COSH}*/
#else /* less precision or not a*/
/* BC3 and not a MC6 ? */
#define FABS(x) (REAL)fabs((double)(x)) /* declare names of basic */
#ifdef LONG_DOUBLE_USED /* floating point */
#define SQRT(x) sqrtlong((x)) /* functions that can be */
/*.IX{SQRT}*/
#else /* used in each of the */
#define SQRT(x) (REAL)sqrt((double)(x)) /* three precisions */
#endif
#define POW(x, y) (REAL)pow((double)(x), \
/*.IX{POW}*/ \
(double)(y))
#define SIN(x) (REAL)sin((double)(x))
/*.IX{SIN}*/
#define COS(x) (REAL)cos((double)(x))
/*.IX{COS}*/
#define EXP(x) (REAL)exp((double)(x))
/*.IX{EXP}*/
#define LOG(x) (REAL)log((double)(x))
/*.IX{LOG}*/
#define ATAN(x) (REAL)atan((double)(x))
/*.IX{ATAN}*/
#define ACOS(x) (REAL)acos((double)(x))
/*.IX{ACOS}*/
#define COSH(x) (REAL)cosh((double)(x))
/*.IX{COSH}*/
#endif
#undef sign
#undef min
#undef max
#define sign(x, y) (((y) < ZERO) ? -FABS(x) : /* |x| times */ \
/*.IX{sign}*/ \
FABS(x)) /* sign of y */
#define min(a, b) (((a) < (b)) ? (a) : (b))
/*.IX{min}*/
#define max(a, b) (((a) > (b)) ? (a) : (b))
/*.IX{max}*/
#define SWAP(typ, a, b) /* swap two objects of */ \
/*.IX{SWAP}*/ \
{ typ temp; temp = a; a = b; b = temp; } /* arbitrary type */
/* ------------------ Macros by Albert Becker ----------------------- */
#define ABS(X) (((X) >= ZERO) ? (X) : -(X)) /* Absolute value of X */
/*.IX{ABS}*/
#define SIGN(X,Y) \
/*.IX{SIGN}*/ \
(((Y) < ZERO) ? -ABS(X) : ABS(X)) /* sign of Y times */
/* ABS(X) */
#define SQR(X) ((X) * (X)) /* square of X */
/*.IX{SQR}*/
#define FORMAT_IN "%lg" /* Input format for REAL */
/*.IX{FORMAT\unt IN}*/
#define FORMAT_LF "% "LZP"f " /* Format l for REAL */
/*.IX{FORMAT\unt LF}*/
#define FORMAT_126LF "% 12.6"LZP"f " /* Format 12.6f for REAL */
/*.IX{FORMAT\unt 126LF}*/
#define FORMAT_2010LF "% 20.10"LZP"f " /* Format 20.10f for REAL */
/*.IX{FORMAT\unt 2010LF}*/
#define FORMAT_2016LF "% 20.16"LZP"f " /* Format 20.16f for REAL */
/*.IX{FORMAT\unt 2016LF}*/
#define FORMAT_LE "% "LZP"e " /* Format e for REAL */
/*.IX{FORMAT\unt LE}*/
#define FORMAT_2016LE "% 20.16"LZP"e " /* Format 20.16e for REAL */
/*.IX{FORMAT\unt 2016LE}*/
/***********************************************************************
* declare all external functions defined in basis.c *
***********************************************************************/
int basis(void); /* find basis for number representation */
REAL mach_eps(void); /* find machine constant */
REAL epsroot(void); /* find square root of machine constant */
REAL epsquad(void); /* find square of machine constant */
REAL maxroot(void); /* Root of the largest representable number ....*/
REAL posmin(void); /* find smallst positive floating point */
/* number */
REAL pi(void); /* find pi */
REAL exp_1(void); /* find e */
REAL sqr(REAL x); /* square a floating point number */
void fehler_melden /* Write error messages to stdout and stderr ....*/
(
char text[], /* error description ........*/
int fehlernummer, /* Number of error ..........*/
char dateiname[], /* file with error .........*/
int zeilennummer /* file name, row number ....*/
);
int umleiten /* Perhaps redirect stdin or stdout to a file */
(
int argc, /* number of arguments in command line ..*/
char *argv[] /* Vector of arguments ..................*/
); /* error code ...........................*/
void readln(void); /* Skip the remainder of line in stdin */
void getline /* Read one line from stdin ....................*/
(
char kette[], /* Vector with the read text ...........*/
int limit /* maximal length of kette .............*/
);
int intervall /* Find the number for a value inside a partition ...*/
(
int n, /* lenght of partition ..................*/
REAL xwert, /* number whose interval index is wanted */
REAL x[] /* partition ............................*/
); /* Index for xwert ......................*/
REAL horner /* Horner scheme for polynomial evaluations .......*/
(
int n, /* Polynomial degree ......*/
REAL a[], /* Polynomial coefficients */
REAL x /* place of evaluation ....*/
); /* Polynomial value at x ..*/
REAL norm_max /* Find the maximum norm of a REAL vector .........*/
(
REAL vektor[], /* vector .................*/
int n /* length of vector .......*/
); /* Maximum norm ...........*/
REAL skalprod /* standard scalar product of two REAL vectors*/
(
REAL v[], /* 1st vector .....................*/
REAL w[], /* 2nd vector .....................*/
int n /* vector length...................*/
); /* scalar product .................*/
void copy_vector /* copy a REAL vector ........................*/
(
REAL ziel[], /* copied vector ............*/
REAL quelle[], /* original vector ..........*/
int n /* length of vector .........*/
);
/*--------------------------------------------------------------------*
* Basic functions chapter 1 (by Albert Becker) ......................*
*--------------------------------------------------------------------*/
long double sqrtlong (long double x);
int comdiv /* Complex division ..........................*/
(
REAL ar, /* Real part of numerator ..........*/
REAL ai, /* Imaginary part of numerator .....*/
REAL br, /* Real part of denominator ........*/
REAL bi, /* Imaginary part of denominator ...*/
REAL * cr, /* Real part of quotient ...........*/
REAL * ci /* Imaginary part of quotient ......*/
);
REAL comabs /* Complex absolute value ....................*/
(
REAL ar, /* Real part .......................*/
REAL ai /* Imaginary part ..................*/
);
void quadsolv /* Complex quadratic equation ................*/
(
REAL ar, /* second degree coefficient .......*/
REAL ai,
REAL br, /* linear coefficient ..............*/
REAL bi,
REAL cr, /* polynomial constant .............*/
REAL ci,
REAL * tr, /* solution ........................*/
REAL * ti
);
void SetVec /* initialize vector .........................*/
(int n, REAL x[], REAL val);
void CopyVec /* copy vector ...............................*/
(int n, REAL source[], REAL dest[]);
int ReadVec /* read vector from input index starts at zero*/
(FILE *fp, int n, REAL x[]);
int ReadVec1 /* read vector from input index starts at one */
(FILE *fp, int n, REAL x[]);
int WriteVec /* write vector to output index starts at zero*/
(FILE *fp, int n, REAL x[]);
int WriteVec1 /* write vector to output index starts at one */
(FILE *fp, int n, REAL x[]);
void SetMat /* initialize matrix .........................*/
(int m, int n, REAL * a[], REAL val);
void CopyMat /* copy matrix ...............................*/
(int m, int n, REAL * source[], REAL * dest[]);
int ReadMat /* read matrix from input index starts at zero*/
(FILE *fp, int m, int n, REAL * a[]);
int ReadMat1 /* read matrix from input index starts at one */
(FILE *fp, int m, int n, REAL * a[]);
int WriteMat /* write matrix to output index starts at zero*/
(FILE *fp, int m, int n, REAL * mat[]);
int WriteMat1 /* write matrix to output index starts at one */
(FILE *fp, int m, int n, REAL * mat[]);
int WriteHead (FILE *fp, char *s); /* write header to file ........*/
int WriteEnd (FILE *fp); /* write separator to file .......*/
void LogError /* write error message to stdout .............*/
(char *s, int rc, char *file, int line);
#endif
/* -------------------------- END basis.h --------------------------- */
vmbloch.h
// --------------------- DECLARATIONS vmblock.h ------------------------
// Header file of vmblock.cpp
#ifndef VMBLOCK_H_INCLUDED
#define VMBLOCK_H_INCLUDED
/***********************************************************************
* symbolic names that let the user choose the kind of dynamical data *
* structure to be allocated upon calling vmalloc() *
***********************************************************************/
#define VEKTOR 0 /* for a REAL vector */
#define VVEKTOR 1 /* for a vector with elements */
/* of given size */
#define MATRIX 2 /* for a REAL matrix */
#define IMATRIX 3 /* for an int matrix */
#define MMATRIX 4 /* for a matrix of 4x4 matrices */
/* (with elements of type `mat4x4') */
#define UMATRIX 5 /* for a lower triangular matrix */
#define PMATRIX 6 /* for a matrix of points in R3 */
/***********************************************************************
* operations on a vector matrix list *
***********************************************************************/
void *vminit /* create an empty vector/matrix list ...........*/
(
void
); /* address of list ...................*/
void *vmalloc /* create a dynamic vector or matrix ............*/
(
void *vmblock, /* address of a vector/matrix list ...*/
int typ, /* kind of vector or matrix ..........*/
size_t zeilen, /* length (vector) or number of rows .*/
size_t spalten /* number of columns or element size .*/
); /* address of the created object .....*/
bool vmcomplete /* check vector/matrix list for lack of memory ..*/
(
void *vmblock /* address of list ...................*/
); /* no lack of memory? ................*/
void vmfree /* free the memory for a vektor/matrix list .....*/
(
void *vmblock /* address of list ...................*/
);
#endif
/* ------------------------- END vmblock.h -------------------------- */