| 1 | /////////////////////////////////////////////////////////////////////////////////
|
|---|
| 2 | //
|
|---|
| 3 | // Levenberg - Marquardt non-linear minimization algorithm
|
|---|
| 4 | // Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
|
|---|
| 5 | // Institute of Computer Science, Foundation for Research & Technology - Hellas
|
|---|
| 6 | // Heraklion, Crete, Greece.
|
|---|
| 7 | //
|
|---|
| 8 | // This program is free software; you can redistribute it and/or modify
|
|---|
| 9 | // it under the terms of the GNU General Public License as published by
|
|---|
| 10 | // the Free Software Foundation; either version 2 of the License, or
|
|---|
| 11 | // (at your option) any later version.
|
|---|
| 12 | //
|
|---|
| 13 | // This program is distributed in the hope that it will be useful,
|
|---|
| 14 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 15 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 16 | // GNU General Public License for more details.
|
|---|
| 17 | //
|
|---|
| 18 | /////////////////////////////////////////////////////////////////////////////////
|
|---|
| 19 |
|
|---|
| 20 | #ifndef _MISC_H_
|
|---|
| 21 | #define _MISC_H_
|
|---|
| 22 |
|
|---|
| 23 | /* common suffix for LAPACK subroutines. Define empty in case of no prefix. */
|
|---|
| 24 | #define LM_LAPACK_SUFFIX _
|
|---|
| 25 | //#define LM_LAPACK_SUFFIX // define empty
|
|---|
| 26 |
|
|---|
| 27 | /* common prefix for BLAS subroutines. Leave undefined in case of no prefix.
|
|---|
| 28 | * You might also need to modify LM_BLAS_PREFIX below
|
|---|
| 29 | */
|
|---|
| 30 | /* f2c'd BLAS */
|
|---|
| 31 | //#define LM_BLAS_PREFIX f2c_
|
|---|
| 32 | /* C BLAS */
|
|---|
| 33 | //#define LM_BLAS_PREFIX cblas_
|
|---|
| 34 |
|
|---|
| 35 | /* common suffix for BLAS subroutines */
|
|---|
| 36 | //#define LM_BLAS_SUFFIX // define empty if a f2c_ or cblas_ prefix was defined for LM_BLAS_PREFIX above
|
|---|
| 37 | #define LM_BLAS_SUFFIX _ // use this in case of no BLAS prefix
|
|---|
| 38 |
|
|---|
| 39 |
|
|---|
| 40 | #define LCAT_(a, b) #a b
|
|---|
| 41 | #define LCAT(a, b) LCAT_(a, b) // force substitution
|
|---|
| 42 | #define RCAT_(a, b) a #b
|
|---|
| 43 | #define RCAT(a, b) RCAT_(a, b) // force substitution
|
|---|
| 44 |
|
|---|
| 45 | #define LM_MK_LAPACK_NAME(s) LM_ADD_PREFIX(LM_CAT_(s, LM_LAPACK_SUFFIX))
|
|---|
| 46 |
|
|---|
| 47 | #ifdef LM_BLAS_PREFIX
|
|---|
| 48 | #define LM_MK_BLAS_NAME(s) LM_CAT_(LM_BLAS_PREFIX, LM_ADD_PREFIX(LM_CAT_(s, LM_BLAS_SUFFIX)))
|
|---|
| 49 | #else
|
|---|
| 50 | #define LM_MK_BLAS_NAME(s) LM_ADD_PREFIX(LM_CAT_(s, LM_BLAS_SUFFIX))
|
|---|
| 51 | #endif
|
|---|
| 52 |
|
|---|
| 53 |
|
|---|
| 54 | #define __BLOCKSZ__ 32 /* block size for cache-friendly matrix-matrix multiply. It should be
|
|---|
| 55 | * such that __BLOCKSZ__^2*sizeof(LM_REAL) is smaller than the CPU (L1)
|
|---|
| 56 | * data cache size. Notice that a value of 32 when LM_REAL=double assumes
|
|---|
| 57 | * an 8Kb L1 data cache (32*32*8=8K). This is a concervative choice since
|
|---|
| 58 | * newer Pentium 4s have a L1 data cache of size 16K, capable of holding
|
|---|
| 59 | * up to 45x45 double blocks.
|
|---|
| 60 | */
|
|---|
| 61 | #define __BLOCKSZ__SQ (__BLOCKSZ__)*(__BLOCKSZ__)
|
|---|
| 62 |
|
|---|
| 63 | /* add a prefix in front of a token */
|
|---|
| 64 | #define LM_CAT__(a, b) a ## b
|
|---|
| 65 | #define LM_CAT_(a, b) LM_CAT__(a, b) // force substitution
|
|---|
| 66 | #define LM_ADD_PREFIX(s) LM_CAT_(LM_PREFIX, s)
|
|---|
| 67 |
|
|---|
| 68 | #define FABS(x) (((x)>=0.0)? (x) : -(x))
|
|---|
| 69 |
|
|---|
| 70 | #ifdef __cplusplus
|
|---|
| 71 | extern "C" {
|
|---|
| 72 | #endif
|
|---|
| 73 |
|
|---|
| 74 | /* blocking-based matrix multiply */
|
|---|
| 75 | extern void slevmar_trans_mat_mat_mult(float *a, float *b, int n, int m);
|
|---|
| 76 | extern void dlevmar_trans_mat_mat_mult(double *a, double *b, int n, int m);
|
|---|
| 77 |
|
|---|
| 78 | /* forward finite differences */
|
|---|
| 79 | extern void slevmar_fdif_forw_jac_approx(void (*func)(float *p, float *hx, int m, int n, void *adata),
|
|---|
| 80 | float *p, float *hx, float *hxx, float delta,
|
|---|
| 81 | float *jac, int m, int n, void *adata);
|
|---|
| 82 | extern void dlevmar_fdif_forw_jac_approx(void (*func)(double *p, double *hx, int m, int n, void *adata),
|
|---|
| 83 | double *p, double *hx, double *hxx, double delta,
|
|---|
| 84 | double *jac, int m, int n, void *adata);
|
|---|
| 85 |
|
|---|
| 86 | /* central finite differences */
|
|---|
| 87 | extern void slevmar_fdif_cent_jac_approx(void (*func)(float *p, float *hx, int m, int n, void *adata),
|
|---|
| 88 | float *p, float *hxm, float *hxp, float delta,
|
|---|
| 89 | float *jac, int m, int n, void *adata);
|
|---|
| 90 | extern void dlevmar_fdif_cent_jac_approx(void (*func)(double *p, double *hx, int m, int n, void *adata),
|
|---|
| 91 | double *p, double *hxm, double *hxp, double delta,
|
|---|
| 92 | double *jac, int m, int n, void *adata);
|
|---|
| 93 |
|
|---|
| 94 | /* e=x-y and ||e|| */
|
|---|
| 95 | extern float slevmar_L2nrmxmy(float *e, float *x, float *y, int n);
|
|---|
| 96 | extern double dlevmar_L2nrmxmy(double *e, double *x, double *y, int n);
|
|---|
| 97 |
|
|---|
| 98 | /* covariance of LS fit */
|
|---|
| 99 | extern int slevmar_covar(float *JtJ, float *C, float sumsq, int m, int n);
|
|---|
| 100 | extern int dlevmar_covar(double *JtJ, double *C, double sumsq, int m, int n);
|
|---|
| 101 |
|
|---|
| 102 | /* box constraints consistency check */
|
|---|
| 103 | extern int slevmar_box_check(float *lb, float *ub, int m);
|
|---|
| 104 | extern int dlevmar_box_check(double *lb, double *ub, int m);
|
|---|
| 105 |
|
|---|
| 106 | /* Cholesky */
|
|---|
| 107 | extern int slevmar_chol(float *C, float *W, int m);
|
|---|
| 108 | extern int dlevmar_chol(double *C, double *W, int m);
|
|---|
| 109 |
|
|---|
| 110 | #ifdef __cplusplus
|
|---|
| 111 | }
|
|---|
| 112 | #endif
|
|---|
| 113 |
|
|---|
| 114 | #endif /* _MISC_H_ */
|
|---|
