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

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

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

revision 1033 by johnpye, Mon Nov 6 07:49:06 2006 UTC revision 1034 by johnpye, Thu Jan 4 05:37:55 2007 UTC
# Line 530  static int savlinnum=0; Line 530  static int savlinnum=0;
530      Evaluate the objective function.      Evaluate the objective function.
531  */  */
532  static boolean calc_objective( slv3_system_t sys){  static boolean calc_objective( slv3_system_t sys){
533    calc_ok = TRUE;    boolean calc_ok = TRUE;
534    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
535    sys->objective = (sys->obj ? relman_eval(sys->obj,&calc_ok,SAFE_CALC) : 0.0);    sys->objective = (sys->obj ? relman_eval(sys->obj,&calc_ok,SAFE_CALC) : 0.0);
536    Asc_SignalHandlerPop(SIGFPE,SIG_IGN);    Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
# Line 548  static boolean calc_objectives( slv3_sys Line 548  static boolean calc_objectives( slv3_sys
548    rfilter.matchvalue =(REL_INCLUDED);    rfilter.matchvalue =(REL_INCLUDED);
549    rlist = slv_get_solvers_obj_list(SERVER);    rlist = slv_get_solvers_obj_list(SERVER);
550    len = slv_get_num_solvers_objs(SERVER);    len = slv_get_num_solvers_objs(SERVER);
551    calc_ok = TRUE;    boolean calc_ok = TRUE;
552    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
553    for (i = 0; i < len; i++) {    for (i = 0; i < len; i++) {
554      if (rel_apply_filter(rlist[i],&rfilter)) {      if (rel_apply_filter(rlist[i],&rfilter)) {
# Line 578  static boolean calc_inequalities( slv3_s Line 578  static boolean calc_inequalities( slv3_s
578    rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);    rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
579    rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);    rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
580    
581    calc_ok = TRUE;    boolean calc_ok = TRUE;
582    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
583    for (rp=sys->rlist;*rp != NULL; rp++) {    for (rp=sys->rlist;*rp != NULL; rp++) {
584      if (rel_apply_filter(*rp,&rfilter)) {      if (rel_apply_filter(*rp,&rfilter)) {
# Line 593  static boolean calc_inequalities( slv3_s Line 593  static boolean calc_inequalities( slv3_s
593    
594  /**  /**
595      Calculates all of the residuals in the current block and computes      Calculates all of the residuals in the current block and computes
596      the residual norm for block status.  Returns true iff calculations      the residual norm for block status.  
597      preceded without error.  
598        @return 0 on failure, non-zero on success
599  */  */
600  static boolean calc_residuals( slv3_system_t sys){  static boolean calc_residuals( slv3_system_t sys){
601    int32 row;    int32 row;
# Line 603  static boolean calc_residuals( slv3_syst Line 604  static boolean calc_residuals( slv3_syst
604    
605    if( sys->residuals.accurate ) return TRUE;    if( sys->residuals.accurate ) return TRUE;
606    
607    calc_ok = TRUE;    boolean calc_ok = TRUE;
608      boolean calc_ok_1;
609    row = sys->residuals.rng->low;    row = sys->residuals.rng->low;
610    time0=tm_cpu_time();    time0=tm_cpu_time();
611    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);    Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
# Line 619  static boolean calc_residuals( slv3_syst Line 621  static boolean calc_residuals( slv3_syst
621        );        );
622      }      }
623  #endif  #endif
624      sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);      sys->residuals.vec[row] = relman_eval(rel,&calc_ok_1,SAFE_CALC);
625        if(!calc_ok_1){
626            calc_ok = FALSE;
627            CONSOLE_DEBUG("error calculating residual for row %d",row);
628        }
629    
630      if (strcmp(CONVOPT,"ABSOLUTE") == 0) {      if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
631        relman_calc_satisfied(rel,FEAS_TOL);        relman_calc_satisfied(rel,FEAS_TOL);
# Line 632  static boolean calc_residuals( slv3_syst Line 638  static boolean calc_residuals( slv3_syst
638    sys->s.block.funcs++;    sys->s.block.funcs++;
639    square_norm( &(sys->residuals) );    square_norm( &(sys->residuals) );
640    sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);    sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
641      if(!calc_ok)CONSOLE_DEBUG("error calculating residuals");
642    return(calc_ok);    return(calc_ok);
643  }  }
644    
# Line 1341  static int calc_pivots(slv3_system_t sys Line 1348  static int calc_pivots(slv3_system_t sys
1348          ERROR_REPORTER_START_HERE(ASC_PROG_ERROR);          ERROR_REPORTER_START_HERE(ASC_PROG_ERROR);
1349          FPRINTF(stderr,"Relation '");          FPRINTF(stderr,"Relation '");
1350          print_rel_name(stderr,sys,rel);          print_rel_name(stderr,sys,rel);
1351          FPRINTF(stderr,"' not pivoted.\n");          FPRINTF(stderr,"' not pivoted.");
1352          error_reporter_end_flush();          error_reporter_end_flush();
1353    
1354          /*          /*
# Line 2565  static void move_to_next_block( slv3_sys Line 2572  static void move_to_next_block( slv3_sys
2572    
2573      sys->residuals.accurate = FALSE;      sys->residuals.accurate = FALSE;
2574      if( !(ok = calc_residuals(sys)) ) {      if( !(ok = calc_residuals(sys)) ) {
2575        ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);        /* error_reporter will have been called somewhere else already */
2576        FPRINTF(MIF(sys),        CONSOLE_DEBUG("Residual calculation errors detected in move_to_next_block.");
         "Residual calculation errors detected in move_to_next_block.\n");  
       error_reporter_end_flush();  
2577      }      }
2578      if( SHOW_LESS_IMPT &&      if( SHOW_LESS_IMPT &&
2579          (sys->s.block.current_size >1 ||          (sys->s.block.current_size >1 ||
# Line 2737  static void reorder_new_block(slv3_syste Line 2742  static void reorder_new_block(slv3_syste
2742      converged (or is -1, to start).      converged (or is -1, to start).
2743  */  */
2744  static void find_next_unconverged_block( slv3_system_t sys){  static void find_next_unconverged_block( slv3_system_t sys){
2745     do {  
2746       do{
2747       move_to_next_block(sys);       move_to_next_block(sys);
2748  #if DEBUG  #if DEBUG
2749       debug_out_var_values(stderr,sys);       debug_out_var_values(stderr,sys);
2750       debug_out_rel_residuals(stderr,sys);       debug_out_rel_residuals(stderr,sys);
2751  #endif  #endif
2752     } while( !sys->s.converged && block_feasible(sys) && !OPTIMIZING(sys) );     }while( !sys->s.converged && block_feasible(sys) && !OPTIMIZING(sys));
2753    
2754     reorder_new_block(sys);     reorder_new_block(sys);
2755  }  }
2756    
# Line 3552  static void slv3_update_linsolqr(slv3_sy Line 3559  static void slv3_update_linsolqr(slv3_sy
3559    linsolqr_set_condition_tolerance(sys->J.sys, PIVOT_TOL);    linsolqr_set_condition_tolerance(sys->J.sys, PIVOT_TOL);
3560  }  }
3561    
3562  static void slv3_presolve(slv_system_t server, SlvClientToken asys){  static int slv3_presolve(slv_system_t server, SlvClientToken asys){
3563    struct var_variable **vp;    struct var_variable **vp;
3564    struct rel_relation **rp;    struct rel_relation **rp;
3565    int32 cap, ind;    int32 cap, ind;
# Line 3566  static void slv3_presolve(slv_system_t s Line 3573  static void slv3_presolve(slv_system_t s
3573      ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);      ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
3574      FPRINTF(stderr,"QRSlv::slv3_presolve: Variable list was never set.");      FPRINTF(stderr,"QRSlv::slv3_presolve: Variable list was never set.");
3575      error_reporter_end_flush();      error_reporter_end_flush();
3576      return;      return 1;
3577    }    }
3578    if( sys->rlist == NULL && sys->obj == NULL ) {    if( sys->rlist == NULL && sys->obj == NULL ) {
3579      ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);      ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
3580      FPRINTF(stderr,"QRSlv::slv3_presolve: Relation list and objective never set.");      FPRINTF(stderr,"QRSlv::slv3_presolve: Relation list and objective never set.");
3581      error_reporter_end_flush();      error_reporter_end_flush();
3582      return;      return 1;
3583    }    }
3584    
3585    if(sys->presolved > 0) { /* system has been presolved before */    if(sys->presolved > 0) { /* system has been presolved before */
# Line 3643  static void slv3_presolve(slv_system_t s Line 3650  static void slv3_presolve(slv_system_t s
3650    update_status(sys);    update_status(sys);
3651    iteration_ends(sys);    iteration_ends(sys);
3652    sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;    sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
3653    
3654      return 0;
3655  }  }
3656    
3657  #ifdef THIS_IS_AN_UNUSED_FUNCTION  #ifdef THIS_IS_AN_UNUSED_FUNCTION
# Line 3677  static boolean slv3_change_basis(slv3_sy Line 3686  static boolean slv3_change_basis(slv3_sy
3686  }  }
3687  #endif /* THIS_IS_AN_UNUSED_FUNCTION */  #endif /* THIS_IS_AN_UNUSED_FUNCTION */
3688    
3689  static void slv3_resolve(slv_system_t server, SlvClientToken asys)  static int slv3_resolve(slv_system_t server, SlvClientToken asys){
 {  
3690    struct var_variable **vp;    struct var_variable **vp;
3691    struct rel_relation **rp;    struct rel_relation **rp;
3692    slv3_system_t sys;    slv3_system_t sys;
# Line 3708  static void slv3_resolve(slv_system_t se Line 3716  static void slv3_resolve(slv_system_t se
3716    sys->objective =  MAXDOUBLE/2000.0;    sys->objective =  MAXDOUBLE/2000.0;
3717    
3718    update_status(sys);    update_status(sys);
3719      return 0;
3720  }  }
3721    
3722    
3723  static void slv3_iterate(slv_system_t server, SlvClientToken asys)  static int slv3_iterate(slv_system_t server, SlvClientToken asys){
 {  
3724    slv3_system_t sys;    slv3_system_t sys;
3725    FILE              *mif;    FILE              *mif;
3726    FILE              *lif;    FILE              *lif;
# Line 3728  static void slv3_iterate(slv_system_t se Line 3736  static void slv3_iterate(slv_system_t se
3736    sys = SLV3(asys);    sys = SLV3(asys);
3737    mif = MIF(sys);    mif = MIF(sys);
3738    lif = LIF(sys);    lif = LIF(sys);
3739    if (server == NULL || sys==NULL) return;    if (server == NULL || sys==NULL) return 1;
3740    if (check_system(SLV3(sys))) return;    if (check_system(SLV3(sys))) return 2;
3741    if( !sys->s.ready_to_solve ) {    if( !sys->s.ready_to_solve ) {
3742      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv: Not ready to solve.");      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv: Not ready to solve.");
3743      return;      return 3;
3744    }    }
3745    
3746    if (sys->s.block.current_block==-1) {    if (sys->s.block.current_block==-1) {
3747      find_next_unconverged_block(sys);      find_next_unconverged_block(sys);
3748        if(!sys->s.calc_ok){
3749          CONSOLE_DEBUG("Calculation errors after find_next_unconverged_block");
3750          return 10;
3751        }
3752      update_status(sys);      update_status(sys);
3753      if( RELNOMSCALE == 1 /*|| (strcmp(SCALEOPT,"RELNOM") == 0) ||      if( RELNOMSCALE == 1 /*|| (strcmp(SCALEOPT,"RELNOM") == 0) ||
3754         (strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0)*/ ){         (strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0)*/ ){
3755        calc_relnoms(sys);        calc_relnoms(sys);
3756      }      }
3757      return;      return 0; /* not sure if this is an error? */
3758    }    }
3759    if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||    if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
3760        LIFDS)) {        LIFDS)) {
# Line 3759  static void slv3_iterate(slv_system_t se Line 3771  static void slv3_iterate(slv_system_t se
3771      sys->s.diverged = 1;      sys->s.diverged = 1;
3772      iteration_ends(sys);      iteration_ends(sys);
3773      update_status(sys);      update_status(sys);
3774      return;      return 4;
3775    }    }
3776  #endif  #endif
3777    /*    /*
# Line 3803  static void slv3_iterate(slv_system_t se Line 3815  static void slv3_iterate(slv_system_t se
3815          iteration_ends(sys);          iteration_ends(sys);
3816          find_next_unconverged_block(sys);          find_next_unconverged_block(sys);
3817          update_status(sys);          update_status(sys);
3818          return;          return 0;
3819        }        }
3820        ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Direct solve found numerically impossible equation given variables solved in previous blocks.");
       FPRINTF(mif,"Direct solve found numerically impossible equation given variables solved in previous blocks.\n");  
       error_reporter_end_flush();  
3821      case -1:      case -1:
3822        sys->s.inconsistent = TRUE;        sys->s.inconsistent = TRUE;
3823    
# Line 3821  static void slv3_iterate(slv_system_t se Line 3831  static void slv3_iterate(slv_system_t se
3831    
3832        iteration_ends(sys);        iteration_ends(sys);
3833        update_status(sys);        update_status(sys);
3834        return;        return 5;
3835      }      }
3836    } /* if fails with a 0, go on to newton a 1x1 */    } /* if fails with a 0, go on to newton a 1x1 */
3837    if( !calc_J(sys) ) {    if( !calc_J(sys) ) {
# Line 3860  static void slv3_iterate(slv_system_t se Line 3870  static void slv3_iterate(slv_system_t se
3870      iteration_ends(sys);      iteration_ends(sys);
3871      find_next_unconverged_block(sys);      find_next_unconverged_block(sys);
3872      update_status(sys);      update_status(sys);
3873      return;      return 0;
3874    }    }
3875    calc_phi(sys);    calc_phi(sys);
3876    
# Line 3880  static void slv3_iterate(slv_system_t se Line 3890  static void slv3_iterate(slv_system_t se
3890      sys->s.diverged = TRUE;      sys->s.diverged = TRUE;
3891      iteration_ends(sys);      iteration_ends(sys);
3892      update_status(sys);      update_status(sys);
3893      return;      return 6;
3894    }    }
3895    
3896    /* CONSOLE_DEBUG("calc_newton..."); */    /* CONSOLE_DEBUG("calc_newton..."); */
# Line 3923  static void slv3_iterate(slv_system_t se Line 3933  static void slv3_iterate(slv_system_t se
3933        sys->s.inconsistent = TRUE;        sys->s.inconsistent = TRUE;
3934        iteration_ends(sys);        iteration_ends(sys);
3935        update_status(sys);        update_status(sys);
3936        return;        return 7;
3937      }      }
3938    
3939  /* end of code by AWW */  /* end of code by AWW */
# Line 4001  static void slv3_iterate(slv_system_t se Line 4011  static void slv3_iterate(slv_system_t se
4011        sys->s.diverged = TRUE;        sys->s.diverged = TRUE;
4012        iteration_ends(sys);        iteration_ends(sys);
4013        update_status(sys);        update_status(sys);
4014        return;        return 8;
4015      }      }
4016    
4017    
# Line 4028  static void slv3_iterate(slv_system_t se Line 4038  static void slv3_iterate(slv_system_t se
4038        sys->s.diverged = TRUE;        sys->s.diverged = TRUE;
4039        iteration_ends(sys);        iteration_ends(sys);
4040        update_status(sys);        update_status(sys);
4041        return;        return 9;
4042      }      }
4043    
4044      /**      /**
# Line 4050  static void slv3_iterate(slv_system_t se Line 4060  static void slv3_iterate(slv_system_t se
4060        sys->s.calc_ok = new_ok;        sys->s.calc_ok = new_ok;
4061        iteration_ends(sys);        iteration_ends(sys);
4062        update_status(sys);        update_status(sys);
4063        return;        return 0;
4064      }      }
4065      if( !new_ok ) {      if( !new_ok ) {
4066        previous = oldphi;        previous = oldphi;
# Line 4105  static void slv3_iterate(slv_system_t se Line 4115  static void slv3_iterate(slv_system_t se
4115      find_next_unconverged_block(sys);      find_next_unconverged_block(sys);
4116    }    }
4117    update_status(sys);    update_status(sys);
4118      return 0;
4119  }  }
4120    
4121    
4122  static void slv3_solve(slv_system_t server, SlvClientToken asys)  static void slv3_solve(slv_system_t server, SlvClientToken asys)
4123  {  {
4124      int err = 0;
4125    slv3_system_t sys;    slv3_system_t sys;
4126    sys = SLV3(asys);    sys = SLV3(asys);
4127    if (server == NULL || sys==NULL) return;    if (server == NULL || sys==NULL) return 1;
4128    if (check_system(sys)) return;    if (check_system(sys)) return 1;
4129    while( sys->s.ready_to_solve ) slv3_iterate(server,sys);    while( sys->s.ready_to_solve ) err = err | slv3_iterate(server,sys);
4130      return err;
4131  }  }
4132    
4133  static mtx_matrix_t slv3_get_jacobian(slv_system_t server, SlvClientToken sys)  static mtx_matrix_t slv3_get_jacobian(slv_system_t server, SlvClientToken sys)

Legend:
Removed from v.1033  
changed lines
  Added in v.1034

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22