/[ascend]/trunk/base/generic/solver/ida.c
ViewVC logotype

Diff of /trunk/base/generic/solver/ida.c

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

revision 947 by johnpye, Sat Nov 25 15:28:56 2006 UTC revision 948 by johnpye, Sun Nov 26 01:36:49 2006 UTC
# Line 582  int integrator_ida_djex(long int Neq, re Line 582  int integrator_ida_djex(long int Neq, re
582      /* print vars */      /* print vars */
583      for(i=0; i < blsys->n_y; ++i){      for(i=0; i < blsys->n_y; ++i){
584          varname = var_make_name(blsys->system, blsys->y[i]);          varname = var_make_name(blsys->system, blsys->y[i]);
585          CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));          CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
586          ASC_FREE(varname);          ASC_FREE(varname);
587      }      }
588    
# Line 590  int integrator_ida_djex(long int Neq, re Line 590  int integrator_ida_djex(long int Neq, re
590      for(i=0; i < blsys->n_y; ++i){      for(i=0; i < blsys->n_y; ++i){
591          if(blsys->ydot[i]){          if(blsys->ydot[i]){
592              varname = var_make_name(blsys->system, blsys->ydot[i]);              varname = var_make_name(blsys->system, blsys->ydot[i]);
593              CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yp,i));              CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
594              ASC_FREE(varname);              ASC_FREE(varname);
595          }else{          }else{
596              varname = var_make_name(blsys->system, blsys->y[i]);              varname = var_make_name(blsys->system, blsys->y[i]);
# Line 616  int integrator_ida_djex(long int Neq, re Line 616  int integrator_ida_djex(long int Neq, re
616              i< enginedata->nrels && relptr != NULL;              i< enginedata->nrels && relptr != NULL;
617              ++i, ++relptr              ++i, ++relptr
618      ){      ){
         relname = rel_make_name(blsys->system, *relptr);  
         CONSOLE_DEBUG("RELATION %d '%s'",i,relname);  
         ASC_FREE(relname);  
   
619          /* get derivatives for this particular relation */          /* get derivatives for this particular relation */
620          status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);          status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
         /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */  
   
         for(j=0;j<count;++j){  
             varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);  
             /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */  
             ASC_FREE(varname);  
         }  
621    
622          if(status){          if(status){
623              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(blsys->system, *relptr);
# Line 637  int integrator_ida_djex(long int Neq, re Line 626  int integrator_ida_djex(long int Neq, re
626              break;              break;
627          }          }
628    
629          for(j=0; j < enginedata->nrels; ++j){          /* output what's going on here ... */
630              DENSE_ELEM(Jac,i,j) = 0;          relname = rel_make_name(blsys->system, *relptr);
631          }          CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
632            fprintf(stderr,"%d: '%s': ",i,relname);
633          for(j=0; j < count; ++j){          ASC_FREE(relname);
634              /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */          for(j=0;j<count;++j){
635              varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);              varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
636              if(varname){              var_yindex = blsys->y_id[variables[j]];
637                  CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);              if(var_yindex >=0){
638                  ASC_FREE(varname);                  fprintf(stderr,"  var[%d]='%s'=y[%d]",variables[j],varname,var_yindex);
639              }else{              }else{
640                  CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);                  fprintf(stderr,"  var[%d]='%s'=ydot[%d]",variables[j],varname,-var_yindex-1);
641              }              }
642                            ASC_FREE(varname);
643              var_yindex = blsys->y_id[variables[j] - 1];          }
644              CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex);          fprintf(stderr,"\n");
645    
646            /* zero the Jacobian row */
647            for(j=0; j < enginedata->nrels; ++j){
648                DENSE_ELEM(Jac,i,j) = 0;
649            }
650        
651            /* insert values into the Jacobian row in appropriate spots */
652            for(j=0; j < count; ++j){          
653                var_yindex = blsys->y_id[variables[j]];
654              if(var_yindex >= 0){              if(var_yindex >= 0){
655                    asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
656                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];                  DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
                 CONSOLE_DEBUG("New value (%d,%d) is %f", i, var_yindex, DENSE_ELEM(Jac,i,var_yindex));  
657              }else{              }else{
658                    asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
659                  DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;                  DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
                 CONSOLE_DEBUG("New value (%d,%d) is %f (deriv)", i, -var_yindex-1, DENSE_ELEM(Jac,i,-var_yindex-1));  
660              }              }
661          }          }
662      }      }
663    
664      CONSOLE_DEBUG("PRINTING JAC");      CONSOLE_DEBUG("PRINTING JAC");
665      for(i=0; i < blsys->n_y; ++i){      fprintf(stderr,"\t");
666        for(j=0; j < blsys->n_y; ++j){
667            if(j)fprintf(stderr,"\t");
668            varname = var_make_name(blsys->system,blsys->y[j]);
669            fprintf(stderr,"%11s",varname);
670            ASC_FREE(varname);
671        }
672        fprintf(stderr,"\n");
673        for(i=0; i < enginedata->nrels; ++i){
674            relname = rel_make_name(blsys->system, enginedata->rellist[i]);
675            fprintf(stderr,"%s\t",relname);
676            ASC_FREE(relname);
677    
678          for(j=0; j < blsys->n_y; ++j){          for(j=0; j < blsys->n_y; ++j){
679              fprintf(stderr,"%8.2e\t",DENSE_ELEM(Jac,i,j));              if(j)fprintf(stderr,"\t");
680                fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
681          }          }
682          fprintf(stderr,"\n");          fprintf(stderr,"\n");
683      }      }
684    
685  #if 1  #if 0
686      double *yval = NV_DATA_S(yy);      double *yval = NV_DATA_S(yy);
687      /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,      /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,
688      I've pulled the jacobians straight from idadenx.c (an IDA example) */      I've pulled the jacobians straight from idadenx.c (an IDA example) */
Legend:
Removed from v.947  
changed lines
  Added in v.948
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22