Marsyas
0.6.0-alpha
|
00001 /********************************************************** 00002 * Reduced and modified version of basis.c by J-P Moreau * 00003 * ------------------------------------------------------- * 00004 * Reference for original basis.c: * 00005 * * 00006 * "Numerical Algorithms with C, by Gisela Engeln-Muellges * 00007 * and Frank Uhlig, Springer-Verlag, 1996". * 00008 **********************************************************/ 00009 #include <marsyas/basis.h> 00010 00011 REAL sqr(REAL x) /* square a floating point number */ 00012 /*********************************************************************** 00013 * Compute the square of a floating point number. * 00014 * * 00015 * global names used: * 00016 * ================== * 00017 * REAL * 00018 ***********************************************************************/ 00019 { 00020 return x * x; 00021 } 00022 00023 REAL norm_max /* Find the maximum norm of a REAL vector .........*/ 00024 ( 00025 REAL vektor[], /* vector .................*/ 00026 int n /* length of vector .......*/ 00027 ) /* Maximum norm ...........*/ 00028 00029 /*********************************************************************** 00030 * Return the maximum norm of a [0..n-1] vector v. * 00031 * * 00032 * global names used: * 00033 * ================== * 00034 * REAL, FABS, ZERO * 00035 ***********************************************************************/ 00036 { 00037 REAL norm, /* local max */ 00038 betrag; /* magnitude of a component */ 00039 00040 for (n--, norm = ZERO; n >= 0; n--, vektor++) 00041 if ((betrag = FABS(*vektor)) > norm) 00042 norm = betrag; 00043 00044 return norm; 00045 } 00046 00047 void copy_vector /* copy a REAL vector ........................*/ 00048 ( 00049 REAL ziel[], /* copied vector ............*/ 00050 REAL quelle[], /* original vector ..........*/ 00051 int n /* length of vector .........*/ 00052 ) 00053 /*********************************************************************** 00054 * copy the n elements of the vector quelle into the vector ziel. * 00055 * * 00056 * global names used: * 00057 * ================== * 00058 * REAL * 00059 ***********************************************************************/ 00060 { 00061 for (n--; n >= 0; n--) 00062 *ziel++ = *quelle++; 00063 } 00064 00065 /*********************************************************************** 00066 * read a vector x of size n from input file fp (index begins at zero) * 00067 * * 00068 * global names used: * 00069 * ================== * 00070 * REAL * 00071 ***********************************************************************/ 00072 int ReadVec (FILE *fp, int n, REAL x[]) 00073 { 00074 int i; 00075 REAL tmp; 00076 00077 for (i = 0; i < n; ++i) 00078 { 00079 if (fscanf (fp,FORMAT_IN, &tmp) <= 0) return (-1); 00080 x[i] = (REAL) tmp; 00081 } 00082 00083 return (0); 00084 } 00085 00086 /*********************************************************************** 00087 * read a vector x of size n from input file fp (index begins at one) * 00088 * * 00089 * global names used: * 00090 * ================== * 00091 * REAL * 00092 ***********************************************************************/ 00093 int ReadVec1 (FILE *fp, int n, REAL x[]) 00094 { 00095 int i; 00096 REAL tmp; 00097 00098 for (i = 1; i < n+1; ++i) 00099 { 00100 if (fscanf (fp,FORMAT_IN, &tmp) <= 0) return (-1); 00101 x[i] = (REAL) tmp; 00102 } 00103 00104 return (0); 00105 } 00106 00107 void SetVec (int n, REAL x[], REAL val) 00108 { 00109 int i; 00110 00111 for (i = 0; i < n; ++i) 00112 x[i] = val; 00113 } 00114 00115 int WriteVec (FILE *fp, int n, REAL x[]) 00116 /*====================================================================* 00117 * * 00118 * Put out vector of length n to output file. Index begins at zero. * 00119 * * 00120 *====================================================================* 00121 * * 00122 * Input parameters: * 00123 * ================ * 00124 * * 00125 * n int n; * 00126 * lenght of vector * 00127 * x REAL x[]; * 00128 * vector * 00129 * * 00130 * Return value : * 00131 * ============= * 00132 * = 0 All ok * 00133 * = -1 Error putting out to stdout * 00134 * * 00135 *====================================================================*/ 00136 { 00137 int i; 00138 00139 for (i = 0; i < n; ++i) 00140 if (fprintf (fp,FORMAT_126LF, x[i]) <= 0) return (-1); 00141 if (fprintf (fp,"\n") <= 0) return (-1); 00142 00143 return 0; 00144 } 00145 00146 int WriteVec1 (FILE *fp, int n, REAL x[]) 00147 /*====================================================================* 00148 * * 00149 * Put out vector of length n to output file. Index begins at one. * 00150 * * 00151 *====================================================================* 00152 * * 00153 * Input parameters: * 00154 * ================ * 00155 * * 00156 * n int n; * 00157 * lenght of vector * 00158 * x REAL x[]; * 00159 * vector * 00160 * * 00161 * Return value : * 00162 * ============= * 00163 * = 0 All ok * 00164 * = -1 Error putting out to stdout * 00165 * * 00166 *====================================================================*/ 00167 { 00168 int i; 00169 00170 for (i = 1; i < n+1; ++i) 00171 if (fprintf (fp,FORMAT_126LF, x[i]) <= 0) return (-1); 00172 if (fprintf (fp,"\n") <= 0) return (-1); 00173 00174 return 0; 00175 } 00176 00177 REAL pi() { 00178 return 4.0*atan(1.0); 00179 } 00180 00181 void CopyMat (int m, int n, REAL * source[], REAL * dest[]) 00182 /*====================================================================* 00183 * * 00184 * Copy the m x n matrix source to the m x n matrix dest. * 00185 * * 00186 *====================================================================* 00187 * * 00188 * Input parameters: * 00189 * ================ * 00190 * m int m; ( m > 0 ) * 00191 * numnber of rows of matrix * 00192 * n int n; ( n > 0 ) * 00193 * number of columns of matrix * 00194 * source REAL * source[]; * 00195 * matrix * 00196 * dest REAL * dest[]; * 00197 * matrix to be copied to * 00198 * * 00199 * Output parameter: * 00200 * ================ * 00201 * dest same as above * 00202 * * 00203 * ATTENTION : WE do not allocate storage for dest here. * 00204 * * 00205 *====================================================================*/ 00206 { 00207 int i, j; 00208 00209 for (i = 0; i < m; ++i) 00210 for (j = 0; j < n; j++) 00211 dest[i][j] = source[i][j]; 00212 } 00213 00214 REAL maxroot(void) /* Root of the largest representable number ....*/ 00215 /*********************************************************************** 00216 * Compute the square root of the largest machine number 2 ^ (MAX_EXP/2)* 00217 * if not already done * 00218 * * 00219 * global names used: * 00220 * ================== * 00221 * REAL, boolean, FALSE, TRUE, SQRT, POSMAX * 00222 ***********************************************************************/ 00223 { 00224 static REAL save_maxroot; 00225 static boolean schon_berechnet = FALSE; 00226 REAL faktor; 00227 unsigned long int n; 00228 00229 if (! schon_berechnet) 00230 { 00231 save_maxroot = ONE; 00232 faktor = TWO; 00233 for (n = MAX_EXP / 2; n > 1; n /= 2, faktor *= faktor) 00234 if (n % 2 != 0) 00235 save_maxroot *= faktor; 00236 save_maxroot *= faktor; 00237 schon_berechnet = TRUE; 00238 } 00239 00240 return save_maxroot; 00241 } 00242 00243 00244 void quadsolv /* Complex quadratic equation ................*/ 00245 ( 00246 REAL ar, /* second degree coefficient .......*/ 00247 REAL ai, 00248 REAL br, /* linear coefficient ..............*/ 00249 REAL bi, 00250 REAL cr, /* polynomial constant .............*/ 00251 REAL ci, 00252 REAL * tr, /* solution ........................*/ 00253 REAL * ti 00254 ) 00255 /*====================================================================* 00256 * * 00257 * Compute the least magnitude solution of the quadratic equation * 00258 * a * t**2 + b * t + c = 0. Here a, b, c and t are complex. * 00259 * 2 * 00260 * Formeula used: t = 2c / (-b +/- sqrt (b - 4ac)). * 00261 * This formula is valid for a=0 . * 00262 * * 00263 *====================================================================* 00264 * * 00265 * Input parameters: * 00266 * ================ * 00267 * ar, ai coefficient of t**2 REAL ar, ai; * 00268 * br, bi coefficient of t REAL br, bi; * 00269 * cr, ci constant term REAL cr, ci; * 00270 * * 00271 * Output parameter: * 00272 * ================ * 00273 * tr, ti complex solution of minimal magnitude * 00274 * REAL *tr, *ti; * 00275 * * 00276 * Macro used : SQRT * 00277 * ============ * 00278 * * 00279 *====================================================================*/ 00280 { 00281 REAL pr, pi, qr, qi, h; 00282 00283 pr = br * br - bi * bi; 00284 pi = TWO * br * bi; /* p = b * b */ 00285 00286 qr = ar * cr - ai * ci; 00287 qi = ar * ci + ai * cr; /* q = a * c */ 00288 00289 pr = pr - (REAL)4.0 * qr; 00290 pi = pi - (REAL)4.0 * qi; /* p = b * b - 4 * a * c */ 00291 00292 h = SQRT (pr * pr + pi * pi); /* q = sqrt (p) */ 00293 00294 qr = h + pr; 00295 if (qr > ZERO) 00296 qr = SQRT (qr * HALF); 00297 else 00298 qr = ZERO; 00299 00300 qi = h - pr; 00301 if (qi > ZERO) 00302 qi = SQRT (qi * HALF); 00303 else 00304 qi = ZERO; 00305 00306 if (pi < ZERO) qi = -qi; 00307 00308 h = qr * br + qi * bi; /* p = -b +/- q, choose sign for large */ 00309 /* magnitude p */ 00310 if (h > ZERO) 00311 { 00312 qr = -qr; 00313 qi = -qi; 00314 } 00315 00316 pr = qr - br; 00317 pi = qi - bi; 00318 h = pr * pr + pi * pi; /* t = (2 * c) / p */ 00319 00320 if (h == ZERO) 00321 { 00322 *tr = ZERO; 00323 *ti = ZERO; 00324 } 00325 else 00326 { 00327 *tr = TWO * (cr * pr + ci * pi) / h; 00328 *ti = TWO * (ci * pr - cr * pi) / h; 00329 } 00330 } //quadsolv 00331 00332 00333 REAL comabs /* Complex absolute value ....................*/ 00334 ( 00335 REAL ar, /* Real part .......................*/ 00336 REAL ai /* Imaginary part ..................*/ 00337 ) 00338 /*====================================================================* 00339 * * 00340 * Complex absolute value of a * 00341 * * 00342 *====================================================================* 00343 * * 00344 * Input parameters: * 00345 * ================ * 00346 * ar,ai REAL ar, ai; * 00347 * Real, imaginary parts of a * 00348 * * 00349 * Return value : * 00350 * ============= * 00351 * Absolute value of a * 00352 * * 00353 * Macros used : SQRT, ABS, SWAP * 00354 * ============= * 00355 * * 00356 *====================================================================*/ 00357 { 00358 if (ar == ZERO && ai == ZERO) return (ZERO); 00359 00360 ar = ABS (ar); 00361 ai = ABS (ai); 00362 00363 if (ai > ar) /* Switch ai and ar */ 00364 SWAP (REAL, ai, ar) 00365 00366 return ((ai == ZERO) ? (ar) : (ar * SQRT (ONE + ai / ar * ai / ar))); 00367 } 00368 00369 00370 int comdiv /* Complex division ..........................*/ 00371 ( 00372 REAL ar, /* Real part of numerator ..........*/ 00373 REAL ai, /* Imaginary part of numerator .....*/ 00374 REAL br, /* Real part of denominator ........*/ 00375 REAL bi, /* Imaginary part of denominator ...*/ 00376 REAL * cr, /* Real part of quotient ...........*/ 00377 REAL * ci /* Imaginary part of quotient ......*/ 00378 ) 00379 /*====================================================================* 00380 * * 00381 * Complex division c = a / b * 00382 * * 00383 *====================================================================* 00384 * * 00385 * Input parameters: * 00386 * ================ * 00387 * ar,ai REAL ar, ai; * 00388 * Real, imaginary parts of numerator * 00389 * br,bi REAL br, bi; * 00390 * Real, imaginary parts of denominator * 00391 * * 00392 * Output parameters: * 00393 * ================== * 00394 * cr,ci REAL *cr, *ci; * 00395 * Real , imaginary parts of the quotient * 00396 * * 00397 * Return value : * 00398 * ============= * 00399 * = 0 ok * 00400 * = 1 division by 0 * 00401 * * 00402 * Macro used : ABS * 00403 * ============ * 00404 * * 00405 *====================================================================*/ 00406 { 00407 REAL tmp; 00408 00409 if (br == ZERO && bi == ZERO) return (1); 00410 00411 if (ABS (br) > ABS (bi)) 00412 { 00413 tmp = bi / br; 00414 br = tmp * bi + br; 00415 *cr = (ar + tmp * ai) / br; 00416 *ci = (ai - tmp * ar) / br; 00417 } 00418 else 00419 { 00420 tmp = br / bi; 00421 bi = tmp * br + bi; 00422 *cr = (tmp * ar + ai) / bi; 00423 *ci = (tmp * ai - ar) / bi; 00424 } 00425 00426 return (0); 00427 } 00428 00429 int ReadMat (FILE *fp, int m, int n, REAL * a[]) 00430 /*====================================================================* 00431 * * 00432 * Read an m x n matrix from input file. Index starts at zero. * 00433 * * 00434 *====================================================================* 00435 * * 00436 * Input parameters : * 00437 * ================== * 00438 * m int m; ( m > 0 ) * 00439 * number of rows of matrix * 00440 * n int n; ( n > 0 ) * 00441 * column number of matrix * 00442 * a REAL * a[]; * 00443 * matrix * 00444 * * 00445 * Output parameter: * 00446 * ================ * 00447 * a matrix * 00448 * * 00449 * ATTENTION : WE do not allocate storage for a here. * 00450 * * 00451 *====================================================================*/ 00452 { 00453 int i, j; 00454 double x; 00455 00456 for (i = 0; i < m; ++i) 00457 for (j = 0; j < n; j++) 00458 { 00459 if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1); 00460 a[i][j] = (REAL) x; 00461 } 00462 00463 return (0); 00464 } 00465 00466 //same as ReadMat, but beginning at index 1 00467 int ReadMat1 (FILE *fp, int m, int n, REAL * a[]) 00468 { 00469 int i, j; 00470 double x; 00471 00472 for (i = 1; i < m+1; ++i) 00473 for (j = 1; j < n+1; j++) 00474 { 00475 if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1); 00476 a[i][j] = (REAL) x; 00477 } 00478 00479 return (0); 00480 } 00481 00482 int WriteMat (FILE *fp, int m, int n, REAL * a[]) 00483 /*====================================================================* 00484 * * 00485 * Put out an m x n matrix in output file. Index starts at zero. * 00486 * * 00487 *====================================================================* 00488 * * 00489 * Input parameters: * 00490 * ================ * 00491 * m int m; ( m > 0 ) * 00492 * row number of matrix * 00493 * n int n; ( n > 0 ) * 00494 * column number of matrix * 00495 * a REAL * a[]; * 00496 * matrix * 00497 * * 00498 * Return value : * 00499 * ============= * 00500 * = 0 all put out * 00501 * = -1 Error writing onto stdout * 00502 * * 00503 *====================================================================*/ 00504 { 00505 int i, j; 00506 00507 if (fprintf (fp,"\n") <= 0) return (-1); 00508 00509 for (i = 0; i < m; ++i) 00510 { 00511 for (j = 0; j < n; j++) 00512 if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1); 00513 00514 if (fprintf (fp,"\n") <= 0) return (-1); 00515 } 00516 if (fprintf (fp,"\n") <= 0) return (-1); 00517 00518 return (0); 00519 } 00520 00521 //same as WriteMat, but beginning at index 1 00522 int WriteMat1 (FILE *fp, int m, int n, REAL * a[]) 00523 { 00524 int i, j; 00525 00526 if (fprintf (fp,"\n") <= 0) return (-1); 00527 00528 for (i = 1; i < m+1; ++i) 00529 { 00530 for (j = 1; j < n+1; j++) 00531 if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1); 00532 00533 if (fprintf (fp,"\n") <= 0) return (-1); 00534 } 00535 if (fprintf (fp,"\n") <= 0) return (-1); 00536 00537 return (0); 00538 } 00539 00540 00541 void SetMat (int m, int n, REAL * a[], REAL val) 00542 /*====================================================================* 00543 * * 00544 * Initialize an m x n matrix with a constant value val . * 00545 * * 00546 *====================================================================* 00547 * * 00548 * Input parameters: * 00549 * ================ * 00550 * m int m; ( m > 0 ) * 00551 * row number of matrix * 00552 * n int n; ( n > 0 ) * 00553 * column number of matrix * 00554 * a REAL * a[]; * 00555 * matrix * 00556 * val constant value * 00557 * * 00558 * Output parameter: * 00559 * ================ * 00560 * a matrix with constant value val in every position * 00561 * * 00562 *====================================================================*/ 00563 { 00564 int i, j; 00565 00566 for (i = 0; i < m; ++i) 00567 for (j = 0; j < n; j++) 00568 a[i][j] = val; 00569 } 00570 00571 static char Separator[] = 00572 "--------------------------------------------------------------------"; 00573 00574 int WriteHead (FILE *fp, char * string) 00575 /*====================================================================* 00576 * * 00577 * Put out header with text in string in file fp. * 00578 * * 00579 *====================================================================* 00580 * * 00581 * Input parameters: * 00582 * ================ * 00583 * string char *string; * 00584 * text of headertext (last byte is a 0) * 00585 * * 00586 * Return value : * 00587 * ============= * 00588 * = 0 All ok * 00589 * = -1 Error writing onto stdout * 00590 * = -2 Invalid text for header * 00591 * * 00592 *====================================================================*/ 00593 { 00594 if (string == NULL) return (-2); 00595 00596 if (fprintf (fp,"\n%s\n%s\n%s\n\n", Separator, string, Separator) <= 0) 00597 return (-1); 00598 00599 return 0; 00600 } 00601 00602 int WriteHead1( char * string) 00603 /*====================================================================* 00604 * * 00605 * Put out header with text in string in stdout. * 00606 * * 00607 *====================================================================* 00608 * * 00609 * Input parameters: * 00610 * ================ * 00611 * string char *string; * 00612 * text of headertext (last byte is a 0) * 00613 * * 00614 * Return value : * 00615 * ============= * 00616 * = 0 All ok * 00617 * = -1 Error writing onto stdout * 00618 * = -2 Invalid text for header * 00619 * * 00620 *====================================================================*/ 00621 { 00622 if (string == NULL) return (-2); 00623 00624 if (printf ("\n%s\n%s\n%s\n\n", Separator, string, Separator) <= 0) 00625 return (-1); 00626 00627 return 0; 00628 } 00629 00630 int WriteEnd (FILE *fp) 00631 /*====================================================================* 00632 * * 00633 * Put out end of writing onto file fp. * 00634 * * 00635 *====================================================================* 00636 * * 00637 * Return value : * 00638 * ============= * 00639 * = 0 All ok * 00640 * = -1 error writing onto stdout * 00641 * * 00642 *====================================================================*/ 00643 { 00644 if (fprintf (fp,"\n%s\n\n", Separator) <= 0) return (-1); 00645 return 0; 00646 } 00647 00648 int WriteEnd1(void) 00649 /*====================================================================* 00650 * * 00651 * Put out end of writing onto stdout. * 00652 * * 00653 *====================================================================* 00654 * * 00655 * Return value : * 00656 * ============= * 00657 * = 0 All ok * 00658 * = -1 error writing onto stdout * 00659 * * 00660 *====================================================================*/ 00661 { 00662 if (printf ("\n%s\n\n", Separator) <= 0) return (-1); 00663 return 0; 00664 } 00665 00666 void LogError (char * string, int rc, char * file, int line) 00667 /*====================================================================* 00668 * * 00669 * Put out error message onto stdout. * 00670 * * 00671 *====================================================================* 00672 * * 00673 * Input parameters: * 00674 * ================ * 00675 * string char *string; * 00676 * text of error massage (final byte is 0) * 00677 * rc int rc; * 00678 * error code * 00679 * file char *file; * 00680 * name of C file in which error was encountered * 00681 * line int line; * 00682 * line number of C file with error * 00683 * * 00684 *====================================================================*/ 00685 { 00686 if (string == NULL) 00687 { 00688 printf ("Unknown ERROR in file %s at line %d\n", file, line); 00689 return; 00690 } 00691 00692 if (rc == 0) 00693 printf ("ERROR: %s, File %s, Line %d\n", string, file, line); 00694 else 00695 printf ("ERROR: %s, rc = %d, File %s, Line %d\n", 00696 string, rc, file, line); 00697 return; 00698 } 00699 00700 00701 REAL epsroot(void) /* Compute square root of the machine constant ...*/ 00702 /*********************************************************************** 00703 * Compute square root of the machine constant, if not already done. * 00704 * * 00705 * global names used: * 00706 * ================== * 00707 * REAL, boolean, FALSE, TRUE, SQRT, MACH_EPS * 00708 ***********************************************************************/ 00709 00710 { 00711 static REAL save_mach_eps_root; 00712 static boolean schon_berechnet = FALSE; 00713 00714 if (! schon_berechnet) 00715 schon_berechnet = TRUE, 00716 save_mach_eps_root = SQRT(MACH_EPS); 00717 00718 return save_mach_eps_root; 00719 } 00720 00721 00722 REAL epsquad(void) /* Find the machine constant squared .........*/ 00723 /*********************************************************************** 00724 * Compute the square of the machine constant, if not already done. * 00725 * * 00726 * global names used: * 00727 * ================== * 00728 * REAL, boolean, FALSE, TRUE, MACH_EPS * 00729 ***********************************************************************/ 00730 00731 { 00732 static REAL save_mach_eps_quad; 00733 static boolean schon_berechnet = FALSE; 00734 00735 if (! schon_berechnet) 00736 schon_berechnet = TRUE, 00737 save_mach_eps_quad = MACH_EPS * MACH_EPS; 00738 00739 return save_mach_eps_quad; 00740 } 00741 00742 // End of file basis_r.cpp