Changeset 9638


Ignore:
Timestamp:
Feb 15, 2017 12:53:58 PM (3 years ago)
Author:
Christian Andersson
Message:

Minor cleanup and removed many direct usage of direct accesses to DlsMat type. Related to ticket:4815

Location:
trunk/RuntimeLibrary/src/jmi
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/RuntimeLibrary/src/jmi/jmi_block_solver.c

    r9569 r9638  
    875875            } else if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) {
    876876                int j, maxInJacIndex = 0;
    877                 realtype **jac = block->J_scale->cols;
    878                 realtype maxInJac = RAbs(jac[0][i]);
     877                realtype maxInJac = RAbs(DENSE_ELEM(block->J_scale, i, 0));
    879878                block->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */
    880879                for(j = 0; j < N; j++) {
    881                     if(RAbs(maxInJac) < RAbs(jac[j][i])) {
    882                         maxInJac = jac[j][i];
     880                    if(RAbs(maxInJac) < RAbs(DENSE_ELEM(block->J_scale, i, j))) {
     881                        maxInJac = DENSE_ELEM(block->J_scale, i, j);
    883882                        maxInJacIndex = j;
    884883                    }
     
    892891            else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) {
    893892                int j, maxInJacIndex = 0;
    894                 realtype **jac = block->J_scale->cols;
    895                 realtype maxInJac = RAbs(jac[0][i]);
     893                realtype maxInJac = RAbs(DENSE_ELEM(block->J_scale, i, 0));
    896894                /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */
    897895                for(j = 0; j < N; j++) {
    898                     if(RAbs(maxInJac) < RAbs(jac[j][i])) {
    899                         maxInJac = jac[j][i];
     896                    if(RAbs(maxInJac) < RAbs(DENSE_ELEM(block->J_scale, i, j))) {
     897                        maxInJac = DENSE_ELEM(block->J_scale, i, j);
    900898                        maxInJacIndex = j;
    901899                    }
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

    r9636 r9638  
    6666/* Calculate matrix vector product vv = M*v or vv = M^T*v */
    6767static void jmi_kinsol_calc_Mv(DlsMat M, int transpose, N_Vector v, N_Vector vv) {
    68     realtype **m = M->cols;
    6968    realtype*  vd = N_VGetArrayPointer(v);
    7069    realtype*  vvd = N_VGetArrayPointer(vv);
    71     long int N = M->N;
     70    long int N = NV_LENGTH_S(v);
    7271    long int i,j;
    7372
     
    7675            vvd[i] = 0;
    7776            for (j=0;j<N;j++){
    78                 vvd[i] += m[i][j] * vd[j];
     77                vvd[i] += DENSE_ELEM(M, j, i) * vd[j];
    7978            }
    8079        }
     
    8584        for (i=0;i<N;i++) {
    8685            for (j=0;j<N;j++){
    87                 vvd[j] += m[i][j] * vd[i];
     86                vvd[j] += DENSE_ELEM(M, j, i) * vd[i];
     87               
    8888            }
    8989        }
     
    687687                (J->data)[i*N+j] = dres;
    688688            }
    689             J->cols[i] = &(J->data)[i*N];
     689            DENSE_COL(J,i) = &(J->data)[i*N];
    690690            block->dx[i] = 0;
    691691        }       
     
    984984            jmi_log_reals(log, topnode, logInfo, "ivs", N_VGetArrayPointer(kin_mem->kin_uu), block->n);
    985985            if(solver->iterationProgressFlag || nniters == 0) {
    986                 jmi_log_fmt(log, topnode, logInfo, "<scaled_residual_norm:%E>", kin_mem->kin_fnorm);
     986                realtype fnorm;
     987                KINGetFuncNorm(solver->kin_mem, &fnorm);
     988               
     989                jmi_log_fmt(log, topnode, logInfo, "<scaled_residual_norm:%E>", fnorm);
    987990                if (block->callbacks->log_options.log_level >= 5) {
    988991                    jmi_log_node_t node;
     
    10031006                    jmi_log_fmt(log, topnode, logInfo, "<max_scaled_residual_index:%I>", max_index);
    10041007                }
    1005                 if (solver->last_fnorm < kin_mem->kin_fnorm && nniters > 0) { /* This is not ment to happen */
     1008                if (solver->last_fnorm < fnorm && nniters > 0) { /* This is not ment to happen */
    10061009                    jmi_log_node_t warning_node_top;
    10071010                    jmi_log_node_t warning_node;
    10081011                    realtype* last_f = N_VGetArrayPointer(solver->last_residual);
    10091012                    warning_node_top = jmi_log_enter_fmt(log, logWarning, "ResidualIncreaseAfterLineSearch", "The residual L2 norm has increased, from <norm_old: %E> to <norm_new: %E>",
    1010                         solver->last_fnorm, kin_mem->kin_fnorm);
     1013                        solver->last_fnorm, fnorm);
    10111014                    warning_node = jmi_log_enter_index_vector_(log, warning_node_top, logWarning,"increased_residuals", 'R');
    10121015                    for (i=0; i<block->n; i++) { /* Go through the residuals and log which ones increased */
     
    10181021                    jmi_log_leave(log, warning_node_top);
    10191022                }
    1020                 solver->last_fnorm = kin_mem->kin_fnorm;
     1023                solver->last_fnorm = fnorm;
    10211024                solver->last_max_residual= max_residual;
    10221025                solver->last_max_residual_index = max_index;
     
    16151618    jmi_kinsol_solver_t* solver = block->solver;
    16161619    int i,j,k;
    1617     realtype **JTJ_c =  solver->JTJ->cols;
    1618     realtype **jac = block->J->cols;
    16191620    int N = block->n;
    16201621    realtype * uscale_data = N_VGetArrayPointer(solver->kin_y_scale);
     
    16231624    for (i=0;i<N;i++) {
    16241625        /* Add the regularization parameter on the diagonal.   */       
    1625         JTJ_c[i][i] = uscale_data[i]*uscale_data[i];
     1626        DENSE_ELEM(solver->JTJ,i,i) = uscale_data[i]*uscale_data[i];
    16261627        /*Calculate value at RTR(i,i) */
    1627         for (k=0;k<N;k++) JTJ_c[i][i] += jac[i][k]*jac[i][k]*fscale_data[k]*fscale_data[k];
     1628        for (k=0;k<N;k++) DENSE_ELEM(solver->JTJ,i,i) += DENSE_ELEM(block->J,k,i)*DENSE_ELEM(block->J,k,i)*fscale_data[k]*fscale_data[k];
    16281629        for (j=i+1;j<N;j++){
    16291630           
    16301631            /*Calculate value at RTR(i,j) */
    1631             JTJ_c[j][i] = 0;
    1632             for (k=0;k<N;k++) JTJ_c[j][i] += jac[j][k]*jac[i][k]*fscale_data[k]*fscale_data[k];
    1633             JTJ_c[i][j] = JTJ_c[j][i];
     1632            DENSE_ELEM(solver->JTJ,i,j) = 0;
     1633            for (k=0;k<N;k++) DENSE_ELEM(solver->JTJ,i,j) += DENSE_ELEM(block->J,k,j)*DENSE_ELEM(block->J,k,i)*fscale_data[k]*fscale_data[k];
     1634            DENSE_ELEM(solver->JTJ,j,i) = DENSE_ELEM(solver->JTJ,i,j);
    16341635        }
    16351636    }
     
    18321833
    18331834    for(j=0; j < N; j++) {
    1834         realtype *jacCol_j = block->J->cols[j];
     1835        realtype *jacCol_j = DENSE_COL(block->J,j);
    18351836        double tempj = Ith(kin_mem->kin_pp, j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j);
    18361837        for(i = 0; i < N; i++) {
     
    18791880
    18801881    /* work_vector2(i) = sum(k=1:n) step(k)*Jac(i,k) */
    1881     for(i = 0; i < N; i++) {
     1882    for(i = 0; i < N; i++) { 
    18821883        double denom = 0;
    18831884        for(j=0; j < N; j++) {
    1884             denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*block->J->cols[j][i]*block->J->cols[j][i];
     1885            denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*DENSE_ELEM(block->J,i,j)*DENSE_ELEM(block->J,i,j);
    18851886        }
    18861887        Ith(solver->work_vector2,i) = denom;
     
    18881889
    18891890    for(j=0; j < N; j++) {
    1890         realtype *jacCol_j = block->J->cols[j];
     1891        realtype *jacCol_j = DENSE_COL(block->J,j);
    18911892        for(i = 0; i < N; i++) {
    18921893            if( Ith(solver->work_vector2,i) >= UNIT_ROUNDOFF) {
     
    19351936   
    19361937    for(j=0; j < N; j++) {
    1937         realtype *jacCol_j = block->J->cols[j];
     1938        realtype *jacCol_j = DENSE_COL(block->J,j);
    19381939        double tempj1 = Ith(kin_mem->kin_pp,j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j)/denom1;
    19391940        double tempj2 = Ith(solver->work_vector2,j)/denom2;
     
    21622163    }
    21632164    else if(solver->J_is_singular_flag) {
    2164         realtype** jac = block->J->cols;
    21652165        if (N == 1) {
    2166             xd[0] = block->nominal[0] * 0.1 *((bd[0] > 0)?1:-1) * ((jac[0][0] > 0)?1:-1);
     2166            xd[0] = block->nominal[0] * 0.1 *((bd[0] > 0)?1:-1) * ((DENSE_ELEM(block->J,0,0) > 0)?1:-1);
    21672167            ret = 0;
    21682168        } else if (solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
     
    21742174            for (i=0;i<N;i++){
    21752175                xd[i] = 0;
    2176                 for (j=0;j<N;j++) xd[i] += jac[i][j]*bd[j]*fscale_data[j]*fscale_data[j];
    2177             }
    2178 
     2176                for (j=0;j<N;j++) xd[i] += DENSE_ELEM(block->J,j,i)*bd[j]*fscale_data[j]*fscale_data[j];
     2177            }
     2178           
    21792179            gnorm = N_VWL2Norm(x, kin_mem->kin_uscale);
    21802180            if((block->callbacks->log_options.log_level >= 5)) {
     
    26382638        char msg[256];
    26392639        realtype* residual_scaling_factors = NV_DATA_S(block->f_scale);
     2640        long int nni = 0, nfe = 0;
     2641       
     2642        KINGetNumFuncEvals(solver->kin_mem, &nfe);
     2643        KINGetNumNonlinSolvIters(solver->kin_mem, &nni);
    26402644
    26412645        /*kin_f(solver->kin_y, solver->work_vector, block); */
     
    26562660        solver->last_fnorm = N_VWL2Norm(kin_mem->kin_fval, block->f_scale);
    26572661        kin_mem->kin_fnorm = solver->last_fnorm;
    2658         sprintf(msg, "nni = %4ld   nfe = %6ld   fnorm = %26.16g", kin_mem->kin_nni, kin_mem->kin_nfe, solver->last_fnorm);
     2662        sprintf(msg, "nni = %4ld   nfe = %6ld   fnorm = %26.16g", nni, nfe, solver->last_fnorm);
    26592663        kin_info("", "KINSolInit", msg, kin_mem->kin_ih_data);
    26602664        jmi_kinsol_print_progress(block, 2);
     
    28732877    jmi_log_t *log = block->log;
    28742878    jmi_log_node_t node;
     2879    long int nniters = 0;
    28752880   
    28762881    flag = block->F(block->problem_data,block->last_accepted_x, NULL, JMI_BLOCK_WRITE_BACK);
     
    29002905    block->force_rescaling        = solver->saved_state->force_rescaling;
    29012906
     2907    KINGetNumNonlinSolvIters(solver->kin_mem, &nniters);
     2908
    29022909    solver->force_new_J_flag = solver->saved_state->force_new_J_flag;
    29032910    solver->handling_of_singular_jacobian_flag = solver->saved_state->handling_of_singular_jacobian_flag;
    29042911    solver->J_is_singular_flag = solver->saved_state->J_is_singular_flag;
    2905     ((struct KINMemRec *)solver->kin_mem)->kin_nnilset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - solver->saved_state->mbset;
     2912    ((struct KINMemRec *)solver->kin_mem)->kin_nnilset = nniters - solver->saved_state->mbset;
    29062913    memcpy(N_VGetArrayPointer(block->f_scale), N_VGetArrayPointer(solver->saved_state->kin_f_scale), block->n*sizeof(realtype));
    29072914    memcpy(N_VGetArrayPointer(solver->kin_y_scale), N_VGetArrayPointer(solver->saved_state->kin_y_scale), block->n*sizeof(realtype));
     
    29092916   
    29102917    if((block->callbacks->log_options.log_level >= 6)) {
    2911         jmi_log_fmt(block->log, node, logInfo, "<mbset: %d> <nni: %d> < nnilset: %d>", solver->saved_state->mbset, &(((struct KINMemRec *)solver->kin_mem)->kin_nni), &(((struct KINMemRec *)solver->kin_mem)->kin_nnilset));
     2918        jmi_log_fmt(block->log, node, logInfo, "<mbset: %d> <nni: %d> < nnilset: %d>", solver->saved_state->mbset, nniters, &(((struct KINMemRec *)solver->kin_mem)->kin_nnilset));
    29122919        jmi_log_ints(block->log, node, logInfo, "handling_singular", &(solver->handling_of_singular_jacobian_flag), 1);
    29132920        jmi_log_ints(block->log, node, logInfo, "is_singular", &(block->force_rescaling), 1);
     
    29322939        jmi_kinsol_solver_t* solver = block->solver;
    29332940        jmi_log_node_t node;
     2941        long int nniters = 0;
    29342942       
    29352943        flag = block->F(block->problem_data,block->last_accepted_x,block->res,JMI_BLOCK_INITIALIZE);
     
    29602968        solver->saved_state->force_new_J_flag      = solver->force_new_J_flag;
    29612969       
     2970        KINGetNumNonlinSolvIters(solver->kin_mem, &nniters);
     2971       
    29622972        solver->saved_state->handling_of_singular_jacobian_flag = solver->handling_of_singular_jacobian_flag;
    29632973        solver->saved_state->J_is_singular_flag = solver->J_is_singular_flag;
    2964         solver->saved_state->mbset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - ((struct KINMemRec *)solver->kin_mem)->kin_nnilset;
     2974        solver->saved_state->mbset = nniters - ((struct KINMemRec *)solver->kin_mem)->kin_nnilset;
    29652975        memcpy(N_VGetArrayPointer(solver->saved_state->kin_f_scale),N_VGetArrayPointer(block->f_scale), block->n*sizeof(realtype));
    29662976        memcpy(N_VGetArrayPointer(solver->saved_state->kin_y_scale),N_VGetArrayPointer(solver->kin_y_scale), block->n*sizeof(realtype));
     
    29682978       
    29692979        if((block->callbacks->log_options.log_level >= 6)) {
    2970             jmi_log_fmt(block->log, node, logInfo, "<mbset: %d> <nni: %d> < nnilset: %d>", solver->saved_state->mbset, &(((struct KINMemRec *)solver->kin_mem)->kin_nni), &(((struct KINMemRec *)solver->kin_mem)->kin_nnilset));
     2980            jmi_log_fmt(block->log, node, logInfo, "<mbset: %d> <nni: %d> < nnilset: %d>", solver->saved_state->mbset, nniters, &(((struct KINMemRec *)solver->kin_mem)->kin_nnilset));
    29712981            jmi_log_ints(block->log, node, logInfo, "handling_singular", &(solver->handling_of_singular_jacobian_flag), 1);
    29722982            jmi_log_ints(block->log, node, logInfo, "is_singular", &(solver->J_is_singular_flag), 1);
Note: See TracChangeset for help on using the changeset viewer.