/[ascend]/trunk/base/generic/integrator/idaanalyse.c
ViewVC logotype

Diff of /trunk/base/generic/integrator/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1248 by johnpye, Sat Jan 27 00:11:34 2007 UTC revision 1249 by johnpye, Sat Jan 27 04:14:58 2007 UTC
# Line 6  Line 6 
6  # include <solver/diffvars.h>  # include <solver/diffvars.h>
7  # include <solver/slv_stdcalls.h>  # include <solver/slv_stdcalls.h>
8  # include <solver/block.h>  # include <solver/block.h>
9  # include <solver/system_impl.h>  # include <solver/diffvars_impl.h>
10  #include <solver/slvDOF.h>  #include <solver/slvDOF.h>
11  #endif  #endif
12    
13  #define ANALYSE_DEBUG  #define ANALYSE_DEBUG
14    
15    #define VARMSG(MSG) \
16        varname = var_make_name(sys->system,v); \
17        CONSOLE_DEBUG(MSG,varname); \
18        ASC_FREE(varname)
19    
20  static int integrator_ida_check_partitioning(IntegratorSystem *sys);  static int integrator_ida_check_partitioning(IntegratorSystem *sys);
21  static int integrator_ida_check_diffindex(IntegratorSystem *sys);  static int integrator_ida_check_diffindex(IntegratorSystem *sys);
22  static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys);  /* static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys); */
23    
24  /**  /**
25      A var can be non-incident. If it *is* non incident and we're going to      A var can be non-incident. If it *is* non incident and we're going to
# Line 54  static int integrator_ida_check_vars(Int Line 59  static int integrator_ida_check_vars(Int
59      SolverDiffVarSequence seq;      SolverDiffVarSequence seq;
60      int vok;      int vok;
61    
 #define VARMSG(MSG) \  
     varname = var_make_name(sys->system,v); \  
     CONSOLE_DEBUG(MSG,varname); \  
     ASC_FREE(varname)  
   
62      CONSOLE_DEBUG("BEFORE CHECKING VARS");      CONSOLE_DEBUG("BEFORE CHECKING VARS");
63      integrator_ida_analyse_debug(sys,stderr);      integrator_ida_analyse_debug(sys,stderr);
64    
# Line 68  static int integrator_ida_check_vars(Int Line 68  static int integrator_ida_check_vars(Int
68      asc_assert(sys->n_y==0);      asc_assert(sys->n_y==0);
69    
70      /* get the the dervative chains from the system */      /* get the the dervative chains from the system */
71      diffvars = sys->system->diffvars;      diffvars = system_get_diffvars(sys->system);
72      if(diffvars==NULL){      if(diffvars==NULL){
73          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
74          return 1;          return 1;
# Line 213  static int integrator_ida_create_lists(I Line 213  static int integrator_ida_create_lists(I
213      asc_assert(sys->ydot==NULL);      asc_assert(sys->ydot==NULL);
214      asc_assert(sys->y_id== NULL);      asc_assert(sys->y_id== NULL);
215    
     CONSOLE_DEBUG("BEFORE MAKING LISTS");  
     integrator_ida_analyse_debug(sys,stderr);  
   
216      /* get the the dervative chains from the system */      /* get the the dervative chains from the system */
217      diffvars = sys->system->diffvars;      diffvars = system_get_diffvars(sys->system);
218      if(diffvars==NULL){      if(diffvars==NULL){
219          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
220          return 1;          return 1;
# Line 271  static int integrator_ida_create_lists(I Line 268  static int integrator_ida_create_lists(I
268          j++;          j++;
269      }      }
270    
     CONSOLE_DEBUG("AFTER MAKING LISTS");  
     integrator_ida_analyse_debug(sys,stderr);  
   
271      return 0;      return 0;
272  }  }
273    
# Line 327  int integrator_ida_analyse(IntegratorSys Line 321  int integrator_ida_analyse(IntegratorSys
321      CONSOLE_DEBUG("Creating lists");      CONSOLE_DEBUG("Creating lists");
322  #endif  #endif
323    
324        CONSOLE_DEBUG("BEFORE MAKING LISTS");
325        integrator_ida_debug(sys,stderr);
326    
327      res = integrator_ida_create_lists(sys);      res = integrator_ida_create_lists(sys);
328      if(res){      if(res){
329          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
# Line 334  int integrator_ida_analyse(IntegratorSys Line 331  int integrator_ida_analyse(IntegratorSys
331      }      }
332    
333  #ifdef ANALYSE_DEBUG  #ifdef ANALYSE_DEBUG
334      CONSOLE_DEBUG("Checking");      CONSOLE_DEBUG("Checking lists");
335  #endif  #endif
336    
337      asc_assert(sys->y);      asc_assert(sys->y);
338      asc_assert(sys->ydot);      asc_assert(sys->ydot);
339      asc_assert(sys->y_id);      asc_assert(sys->y_id);
340    
341      integrator_ida_analyse_debug(sys,stderr);      integrator_ida_debug(sys,stderr);
342    
343      if(integrator_ida_check_diffindex(sys)){      if(integrator_ida_check_diffindex(sys)){
344          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
# Line 388  int integrator_ida_analyse(IntegratorSys Line 385  int integrator_ida_analyse(IntegratorSys
385  #endif  #endif
386    
387      /* get the the dervative chains from the system */      /* get the the dervative chains from the system */
388      diffvars = sys->system->diffvars;      diffvars = system_get_diffvars(sys->system);
389    
390      /* check the indep var is same as was located elsewhere */      /* check the indep var is same as was located elsewhere */
391      if(diffvars->nindep>1){      if(diffvars->nindep>1){
# Line 402  int integrator_ida_analyse(IntegratorSys Line 399  int integrator_ida_analyse(IntegratorSys
399          return 5;          return 5;
400      }      }
401    
402    #ifdef ANALYSE_DEBUG
403        CONSOLE_DEBUG("Collecting observed variables");
404    #endif
405    
406      /* get the observations */      /* get the observations */
407      /** @TODO this should use a 'problem provider' API instead of static data      /** @TODO this should use a 'problem provider' API instead of static data
408      from the system object */      from the system object */
# Line 414  int integrator_ida_analyse(IntegratorSys Line 415  int integrator_ida_analyse(IntegratorSys
415          ASC_FREE(varname);          ASC_FREE(varname);
416      }      }
417    
     /*   - 'y' list as [ya|yd] */  
     /*   - sparsity pattern for dF/dy and dF/dy' */  
     /*   - sparsity pattern for union of above */  
     /*   - block decomposition based on above */  
     /*   - block decomposition results in reordering of y and y' */  
     /*   - boundaries (optional) */  
     /* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");  
     return -1; */  
   
418      return 0;      return 0;
419  }  }
420    
# Line 530  int integrator_ida_block_check(Integrato Line 522  int integrator_ida_block_check(Integrato
522      return res;      return res;
523  }  }
524    
525    /* return 0 on succes */
526    static int check_dups(struct var_variable **list,int n,int allownull){
527        int i,j;
528        struct var_variable *v;
529        for(i=0; i< n; ++i){
530            v=list[i];
531            if(v==NULL){
532                if(allownull)continue;
533                else return 2;
534            }
535            for(j=0; j<i-1;++j){
536                if(list[j]==NULL)continue;
537                if(v==list[j])return 1;
538            }
539        }
540        return 0;
541    }
542    
543    /*
544        We are going to check that
545    
546          - n_y in range (0,n_vars)
547          - no duplicates anywhere in the varlist
548          - no duplicates in non-NULL elements of ydot
549    
550          - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
551          - var_sindex for first n_y vars match their position the the solver's var list
552    
553          - ydot contains n_ydot non-NULL elements, each match the deriv filter.
554          - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
555    
556          - y_id contains n_ydot elements, with int values in range [0,n_y)
557          - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
558    
559          - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
560    */
561  static int integrator_ida_check_diffindex(IntegratorSystem *sys){  static int integrator_ida_check_diffindex(IntegratorSystem *sys){
562      int i, nv, err = 0;      int i, n_vars, n_ydot=0;
563      struct var_variable **vlist;      struct var_variable **list, *v;
564        char *varname;
565      int diffindex;      int diffindex;
566        const char *msg;
567    
568      CONSOLE_DEBUG("Checking diffindex vector");      CONSOLE_DEBUG("Checking diffindex vector");
569    
570      if(sys->y_id == NULL){      if(sys->y_id == NULL || sys->y == NULL || sys->ydot == NULL){
571          CONSOLE_DEBUG("y_id not allocated");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
572            return 1;
573        }
574        
575        list = slv_get_solvers_var_list(sys->system);
576        n_vars = slv_get_num_solvers_vars(sys->system);
577        
578        if(check_dups(list,n_vars,FALSE)){
579            ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
580          return 1;          return 1;
581      }      }
582    
583      vlist = slv_get_solvers_var_list(sys->system);      if(check_dups(sys->ydot,n_vars,TRUE)){
584      nv = slv_get_num_solvers_vars(sys->system);          ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
585            return 1;
586        }
587    
588      for(i=0; i<nv; ++i){      /* check n_y in range */
589          if(var_deriv(vlist[i])){      if(sys->n_y <= 0 || sys->n_y >= n_vars){
590              if(i < sys->n_y){          ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",sys->n_y, n_vars);
591                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv in wrong position");          return 1;
592                  err++;      }
593              }else{  
594                  if(var_sindex(vlist[i])!=i){      /* check first n_y vars */
595                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv var with incorrect sindex");      for(i=0; i < sys->n_y; ++i){
596                      err++;          v = list[i];
597                  }          if(!var_apply_filter(v, &integrator_ida_nonderiv)){
598                  diffindex = sys->y_id[i - sys->n_y];              msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
599                  if(diffindex >= sys->n_y){          }else if(v != sys->y[i]){
600                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff y_id[%d] value too big",i-sys->n_y);              msg = "'%s' not matched in y vector"; goto finish;
601                      err++;          }else if(var_sindex(v) != i){
602                  }              msg = "'%s' has wrong var_sindex"; goto finish;
603                  if(var_sindex(vlist[diffindex])!=diffindex){          }
604                      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff var with incorrect sindex");          /* meets filter, matches in y vector, has correct var_sindex. */
605                      err++;      }
606                  }  
607              }      /* count non-NULLs in ydot */
608          }else{      for(i=0; i < sys->n_y; ++i){
609              if(i!=var_sindex(vlist[i])){          v = sys->ydot[i];
610                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"var_sindex incorrect");          if(v!=NULL){
611                  err++;              if(var_sindex(v) < sys->n_y){
612                    msg = "'%s' has var_sindex < n_y"; goto finish;
613                }else if(!var_apply_filter(v,&integrator_ida_deriv)){
614                    msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
615              }              }
616                /* lies beyond first n_y vars, match deriv filter */
617                n_ydot++;
618          }          }
619      }      }
620      if(err){  
621          CONSOLE_DEBUG("Errors found in diffindex");      /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
622        for(i=sys->n_y; i< sys->n_y + n_ydot; ++i){
623            v = list[i];
624            if(!var_apply_filter(v,&integrator_ida_deriv)){
625                msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
626            }
627        }
628    
629        /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
630        for(i=0; i<n_ydot; ++i){
631            if(sys->y_id[i] < 0 || sys->y_id[i] >= sys->n_y){
632                ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
633                return 1;
634            }
635            v = sys->ydot[sys->y_id[i]];
636            if(var_sindex(v) - sys->n_y != i){
637                msg = "'%s' not index correctly from y_id"; goto finish;
638            }
639        }
640    
641        /* check remaining vars fail both filters */
642        for(i=sys->n_y + n_ydot; i<n_vars; ++i){
643            v = list[i];
644            if(var_apply_filter(v,&integrator_ida_nonderiv)){
645                msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
646            }
647            if(var_apply_filter(v,&integrator_ida_deriv)){
648                CONSOLE_DEBUG("position = %d",i);
649                msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
650            }
651      }      }
652      return err;  
653        return 0;      
654    finish:
655        varname = var_make_name(sys->system,v);
656        ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
657        ASC_FREE(varname);
658        return 1;          
659  }  }
660    
661  /**  /**
Legend:
Removed from v.1248  
changed lines
  Added in v.1249
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22