Changeset 29934 in project


Ignore:
Timestamp:
10/17/13 03:55:16 (8 years ago)
Author:
Ivan Raikov
Message:

nemo: implemented gsl steady state solver in nest backend

Location:
release/4/nemo/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/nemo-nest.scm

    r29932 r29934  
    323323
    324324
    325 (define fsolve-prelude-header
     325(define kinsol-prelude-header
    326326#<<EOF
    327327  /**
     
    510510)
    511511
    512 (define fsolve-prelude
     512(define (fsolve-prelude method)
     513  (case method
     514    ((cvode ida)
    513515#<<EOF
    514516
    515   int fsolve (KINSysFn f, N_Vector fval, void *user_data, std::string name)
     517  int fsolve (KINSysFn f, int N, N_Vector fval, void *user_data, std::string name)
    516518  {
    517       int status, N;
     519      int status;
    518520      N_Vector u0, sc;
    519521      void *kmem;
    520 
    521       N = NV_LENGTH_S(fval);
    522522
    523523      u0 = N_VNew_Serial(N);
     
    551551
    552552EOF
     553  )
     554  (else
     555#<<EOF
     556
     557  int fsolve (int (*fss)(const gsl_vector *, void *user_data, gsl_vector *), int N, gsl_vector *fval, void *user_data, std::string name)
     558  {
     559     
     560      const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid;
     561      gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, N);
     562      gsl_multiroot_function f = {fss, N, user_data};
     563
     564      int status, iter;
     565      gsl_vector *x = gsl_vector_alloc (N);
     566     
     567      for (int i = 0; i < N; i++)
     568      {
     569         gsl_vector_set (x, i, 0.0);
     570      }
     571
     572      gsl_multiroot_fsolver_set (s, &f, x);
     573
     574      iter = 0;
     575      do
     576      {
     577         iter++;
     578         status = gsl_multiroot_fsolver_iterate (s);
     579
     580         if ((status == GSL_EBADFUNC) ||
     581             (status == GSL_ENOPROG))
     582         {
     583            throw GSLSolverFailure(name, status);
     584         }
     585         
     586         status =  gsl_multiroot_test_residual (s->f, 1e-7);
     587      }
     588      while (status == GSL_CONTINUE && iter < 1000);
     589
     590      for (int i = 0; i < N; i++)
     591      {
     592         gsl_vector_set (fval, i, gsl_vector_get (s->x, i));
     593      }
     594
     595      gsl_vector_free (x);
     596      gsl_multiroot_fsolver_free (s);
     597
     598      return 0;
     599  }
     600
     601EOF
     602   
     603))
    553604)
    554 
    555605
    556606
     
    664714           ("#include <gsl/gsl_sf_exp.h>")
    665715           ("#include <gsl/gsl_odeiv2.h>")
     716           ("#include <gsl/gsl_vector.h>")
     717           ("#include <gsl/gsl_multiroots.h>")
    666718           ("#define Ith(v,i)    (v[i])")
    667719           ))
     
    681733        (case method
    682734          ((cvode ida)
    683            (pp indent ,fsolve-prelude-header))))
     735           (pp indent ,kinsol-prelude-header))))
    684736     
    685737 
     
    719771                              int ,(sprintf "~A_steadystate" sysname)
    720772                              "(N_Vector, N_Vector, void*)" #\;) ,nl)
    721            ))
    722           )
     773           )
     774          (else
     775           (pp indent (extern "\"C\""
     776                              int ,(sprintf "~A_steadystate" sysname)
     777                              "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl)
     778           )
     779          ))
    723780
    724781    (pp indent
     
    806863        ((cvode ida)
    807864         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
    808                             "(N_Vector, N_Vector, void*)" #\;) ,nl)
    809          ))
    810       )
     865                            "(N_Vector, N_Vector, void*)" #\;) ,nl))
     866        (else
     867         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
     868                            "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl))
     869        ))
     870         
    811871
    812872  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
     
    25062566          ((cvode ida)
    25072567
    2508            (begin
     2568           (let ((N (length steady-state-index-map))
     2569                 (ssvect (gensym 'ssvect)))
    25092570
    25102571             (pp indent+
    2511                  ("N_Vector ssvect;")
    2512                  ,(expr->string/C++ (sprintf "N_VNew_Serial(~A);" (length steady-state-index-map)) 'ssvect)
    2513                  (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
     2572                 ,(sprintf "N_Vector ~A;" ssvect)
     2573                 ,(expr->string/C++ (sprintf "N_VNew_Serial(~A)" N) ssvect)
     2574                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
    25142575             
    25152576             (for-each
     
    25202581                  (if si
    25212582                        (pp indent+
    2522                             ,(expr->string/C++ (ith 'ssvect si)
     2583                            ,(expr->string/C++ (ith ssvect si)
    25232584                                               (sprintf "y_[~A]" ni))
    2524                             ,(expr->string/C++ (ith 'ssvect si)
     2585                            ,(expr->string/C++ (ith ssvect si)
    25252586                                               n))
    25262587                        )
     
    25282589              rate-eq-defs)
    25292590             
    2530              (pp indent+ ("N_VDestroy_Serial (ssvect);"))
     2591             (pp indent+ (,(sprintf "N_VDestroy_Serial (~A);" ssvect)))
    25312592
    25322593             ))
     2594
     2595          (else
     2596
     2597           (let ((N (length steady-state-index-map))
     2598                 (ssvect (gensym 'ssvect))
     2599                 )
     2600
     2601             (pp indent+
     2602                 ,(sprintf "gsl_vector *~A;" ssvect)
     2603                 ,(expr->string/C++ (sprintf "gsl_vector_alloc (~A)" N) ssvect)
     2604                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
     2605             
     2606             (for-each
     2607              (lambda (def)
     2608                (let* ((n (first def))
     2609                       (si (lookup-def n steady-state-index-map))
     2610                       (ni (lookup-def n state-index-map)))
     2611                  (if si
     2612                        (pp indent+
     2613                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
     2614                                               (sprintf "y_[~A]" ni))
     2615                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
     2616                                               n))
     2617                        )
     2618                  ))
     2619              rate-eq-defs)
     2620             
     2621             (pp indent+ (,(sprintf "gsl_vector_free (~A);" ssvect)))
     2622
     2623             ))
     2624
    25332625          ))
    25342626   
     
    29042996
    29052997 
    2906   (let* ((N (length steady-state-index-map))
     2998  (let* (
     2999         (N (length steady-state-index-map))
    29073000
    29083001         (c-eqs (lookup-def 'c-eqs constraints))
     
    29573050
    29583051         )
    2959    
    2960     (pp indent ,nl
    2961         (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(N_Vector u, N_Vector f, void* pnode)" ,nl)
    2962         ( #\{))
     3052
     3053    (case method
     3054      ((ida cvode)
     3055       (pp indent ,nl
     3056           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(N_Vector u, N_Vector f, void* pnode)" ,nl)
     3057           ( #\{)))
     3058      (else
     3059       (pp indent ,nl
     3060           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(const gsl_vector *u, void *pnode, gsl_vector *f)" ,nl)
     3061           ( #\{)))
     3062      )
    29633063
    29643064    (pp indent+
     
    29973097     const-defs)
    29983098
    2999     (for-each
    3000      (lambda (def)
    3001        (let* ((n   (first def))
    3002               (ni  (lookup-def n steady-state-index-map)))
    3003          (if ni (pp indent+ ,(expr->string/C++ (ith 'u ni) n)))
    3004          ))
    3005      rate-eq-defs)
     3099    (case method
     3100      ((ida cvode)
     3101       (for-each
     3102        (lambda (def)
     3103          (let* ((n   (first def))
     3104                 (ni  (lookup-def n steady-state-index-map)))
     3105            (if ni (pp indent+ ,(expr->string/C++ (ith 'u ni) n)))
     3106            ))
     3107        rate-eq-defs))
     3108      (else
     3109       (for-each
     3110        (lambda (def)
     3111          (let* ((n   (first def))
     3112                 (ni  (lookup-def n steady-state-index-map)))
     3113            (if ni (pp indent+ ,(expr->string/C++ (sprintf "gsl_vector_get (u, ~A)" ni) n)))
     3114            ))
     3115        rate-eq-defs))
     3116      )
     3117       
    30063118
    30073119    (let ((vi (lookup-def 'v state-index-map)))
     
    30163128                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
    30173129              init-order)
    3018    
    3019     (for-each
    3020      (lambda (def)
    3021        (let* ((n   (first def))
    3022               (ni  (lookup-def n steady-state-index-map))
    3023               (b   (second def))
    3024               (lb  (find-locals (list b))))
    3025          (if ni
    3026              (begin
    3027                (if (not (null? lb))
    3028                    (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
    3029                (pp indent+
    3030                    (,(expr->string/C++ b (ith 'f ni))))))
    3031          ))
    3032      rate-eq-defs)
     3130
     3131    (case method
     3132      ((ida cvode)
     3133       (for-each
     3134        (lambda (def)
     3135          (let* ((n   (first def))
     3136                 (ni  (lookup-def n steady-state-index-map))
     3137                 (b   (second def))
     3138                 (lb  (find-locals (list b))))
     3139            (if ni
     3140                (begin
     3141                  (if (not (null? lb))
     3142                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
     3143                  (pp indent+
     3144                      (,(expr->string/C++ b (ith 'f ni))))))
     3145            ))
     3146        rate-eq-defs))
     3147      (else
     3148       (for-each
     3149        (lambda (def)
     3150          (let* ((n   (first def))
     3151                 (ni  (lookup-def n steady-state-index-map))
     3152                 (b   (second def))
     3153                 (lb  (find-locals (list b))))
     3154            (if ni
     3155                (let ((tmp (gensym 't)))
     3156                  (pp indent+ (double ,tmp #\;))
     3157                  (if (not (null? lb))
     3158                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
     3159                  (pp indent+
     3160                      ,(expr->string/C++ b tmp)
     3161                      (,(sprintf "gsl_vector_set (f,~A,~A);" ni tmp)))
     3162                  ))
     3163            ))
     3164        rate-eq-defs))
     3165      )
    30333166
    30343167    (pp indent ,nl
     
    32323365
    32333366          (if (not (null? steady-state-index-map))
    3234               (case method
    3235                 ((cvode ida)
    3236                  (with-output-to-port cpp-output
    3237                    (lambda ()
    3238                      (pp indent ,fsolve-prelude))))))
     3367              (with-output-to-port cpp-output
     3368                (lambda ()
     3369                  (pp indent ,(fsolve-prelude method)))
     3370                ))
    32393371
    32403372          ;; user-defined functions
  • release/4/nemo/trunk/nemo-version.scm

    r29930 r29934  
    11
    2 (define nemo-version "8.45")
     2(define nemo-version "8.46")
Note: See TracChangeset for help on using the changeset viewer.