Ticket #3879: linear_solution_check_v2.patch

File linear_solution_check_v2.patch, 7.6 KB (added by Christian Andersson, 5 years ago)
  • RuntimeLibrary/src/jmi/jmi_linear_solver.h

     
    3333extern void dgels_(char* TRANS, int* M, int* N, int* NRHS,double* A,int* LDA, double* B,int* LDB,double* WORK,int* LWORK,int* INFO );
    3434extern int dgeequ_(int *m, int *n, double *a, int * lda, double *r__, double *c__, double *rowcnd, double
    3535    *colcnd, double *amax, int *info);
     36extern void dgesdd_(char* JOBZ, int* M, int* N, double* A, int* LDA,
     37            double* S, double* U, int* LDU, double* VT, int* LDVT,
     38            double* WORK, int* LWORK, int* IWORK, int* INFO);
    3639
    3740extern int dlaqge_(int *m, int *n, double *a, int * lda, double *r__, double *c__, double *rowcnd, double
    3841    *colcnd, double *amax, char *equed);
     42extern void dgecon_(char *norm, int *n, double *a, int *lda, double *anorm, double *rcond,
     43             double *work, int *iwork, int *info);
     44extern double dlange_(char *norm, int *m, int *n, double *a, int *lda,
     45             double *work);
     46extern double dlamch_(char *cmach);
    3947
    4048typedef struct jmi_linear_solver_t jmi_linear_solver_t;
    4149
     
    5462    jmi_real_t* factorization;      /**< \brief Matrix for storing the Jacobian factorization */
    5563    jmi_real_t* jacobian;         /**< \brief Matrix for storing the Jacobian */
    5664    jmi_real_t* jacobian_temp;         /**< \brief Matrix for storing the Jacobian */
     65    jmi_real_t* jacobian_temp_augmented;         /**< \brief Matrix for storing the temporary augmented Jacobian */
    5766    jmi_real_t* singular_values;  /**< \brief Vector for the singular values of the Jacobian */
    5867    double* rScale;               /**< \brief Row scaling of the Jacobian matrix */
    5968    double* cScale;               /**< \brief Column scaling of the Jacobian matrix */
     
    6372    int iwork;
    6473    double* zero_vector;
    6574    jmi_real_t* rwork;
     75    int* iwork_array;
     76    double* dgesdd_work;          /**< \brief Work vector for desdd */
     77    int dgesdd_lwork;               /**< \brief Work vector for desdd */
     78    int* dgesdd_iwork;              /**< \brief Work vector for desdd */
    6679};
    6780
    6881#endif /* _JMI_LINEAR_SOLVER_H */
  • RuntimeLibrary/src/jmi/jmi_linear_solver.c

     
    2323#include "jmi_block_solver_impl.h"
    2424#include "jmi_log.h"
    2525
     26
     27/* Estimate condition number utilizing dgecon from LAPACK*/
     28/*
     29static double jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block) {
     30    jmi_linear_solver_t* solver = block->solver;
     31    char norm = 'I';
     32    int N = block->n;
     33    double jnorm = 1.0;
     34    double rcond = 1.0;
     35    int info;
     36
     37    Compute infinity norm of J to be used with dgecon
     38    jnorm = dlange_(&norm, &N, &N, solver->jacobian, &N, solver->rwork);
     39
     40    Compute reciprocal condition number
     41    dgecon_(&norm, &N, solver->factorization, &N, &jnorm, &rcond, solver->rwork, solver->iwork_array,&info);
     42
     43    return 1.0/rcond;
     44}
     45*/
     46
     47
    2648int jmi_linear_solver_new(jmi_linear_solver_t** solver_ptr, jmi_block_solver_t* block) {
    2749    jmi_linear_solver_t* solver= (jmi_linear_solver_t*)calloc(1,sizeof(jmi_linear_solver_t));
    2850/*    jmi_t* jmi = block->jmi;
     
    3860    solver->factorization = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    3961    solver->jacobian = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    4062    solver->jacobian_temp = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
     63    solver->jacobian_temp_augmented = (jmi_real_t*)calloc(n_x*(n_x+1),sizeof(jmi_real_t));
    4164    solver->singular_values = (jmi_real_t*)calloc(n_x,sizeof(jmi_real_t));
    4265    solver->rScale = (double*)calloc(n_x,sizeof(double));
    4366    solver->cScale = (double*)calloc(n_x,sizeof(double));
     
    4669   
    4770    solver->singular_jacobian = 0;
    4871    solver->iwork = 5*n_x;
     72    solver->iwork_array = (int*)calloc(solver->iwork,sizeof(int));
    4973    solver->rwork = (double*)calloc(solver->iwork,sizeof(double));
    5074    solver->zero_vector = (double*)calloc(n_x,sizeof(double));
     75   
     76    solver->dgesdd_lwork = 10*n_x;
     77    solver->dgesdd_work  = (double*)calloc(solver->dgesdd_lwork,sizeof(double));
     78    solver->dgesdd_iwork = (int*)calloc(8*n_x,sizeof(int));
    5179
    5280    for (i=0; i<n_x; i++) {
    5381        solver->zero_vector[i] = 0.0;
     
    164192    }
    165193   
    166194    if((block->callbacks->log_options.log_level >= 5)) {
     195        /* double rcond = jmi_calculate_jacobian_condition_number(block); */
    167196        jmi_log_reals(block->log, destnode, logInfo, "initial_guess", block->initial, block->n);
    168         jmi_log_reals(block->log, destnode, logInfo, "b", block->res, block->n);     
     197        jmi_log_reals(block->log, destnode, logInfo, "b", block->res, block->n);
     198        /* jmi_log_reals(block->log, destnode, logInfo, "condition_number", &rcond, 1); */
    169199        jmi_log_leave(block->log, destnode);
    170200    }
    171201 
     
    174204    i = 1; /* One rhs to solve for */
    175205     
    176206    if (solver->singular_jacobian == 1){
     207        char jobz = 'N', ceps = 'P', csfmin = 'S';
     208        int j,rank_aug = 0, m_x = n_x + 1;
     209        double eps = dlamch_(&ceps);
     210        double sfmin = dlamch_(&csfmin);
     211        double threshold;
     212       
     213        memcpy(solver->jacobian_temp_augmented, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
     214        memcpy(&solver->jacobian_temp_augmented[n_x*n_x], block->res, n_x*sizeof(jmi_real_t));
    177215        /*
    178216         *   DGELSS - compute the minimum norm solution to      a real
    179217         *   linear least squares problem
     
    191229            return -1;
    192230        }
    193231       
     232        /* Verify solution by computing the rank of the augmented Jacobian [A,b] */
     233        /*
     234         *       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
     235         *                          LWORK, IWORK, INFO )
     236         */
     237        dgesdd_(&jobz, &n_x, &m_x, solver->jacobian_temp_augmented, &n_x,
     238                solver->singular_values, NULL, &n_x, NULL, &n_x,
     239                solver->dgesdd_work, &solver->dgesdd_lwork, solver->dgesdd_iwork, &info);
     240               
     241        if(info != 0) {
     242            jmi_log_node(block->log, logError, "Error", "DGESDD failed to compute the SVD of the augmented system [A,b] in <block: %d> with error code <error: %s>", block->label, info);
     243            return -1;
     244        }
     245       
     246        threshold = (eps*solver->singular_values[0] > sfmin)? eps*solver->singular_values[0]: sfmin;
     247        for (j = 0; j < n_x; j++)
     248        {
     249            if (solver->singular_values[j] > threshold) {
     250                rank_aug += 1;
     251            } else {
     252                break;
     253            }
     254        }
     255       
     256        if (rank_aug != rank) {
     257            jmi_log_node(block->log, logError, "Error", "Rank test failed, no solution exists in <block: %d>. Jacobian <rank: %d>, Augmented Jacobian <rank: %d>", block->label, rank, rank_aug);
     258            return -1;
     259        }
     260       
    194261        if(block->callbacks->log_options.log_level >= 5){
    195262            jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
    196263        }
     
    236303    free(solver->jacobian);
    237304    free(solver->singular_values);
    238305    free(solver->jacobian_temp);
     306    free(solver->jacobian_temp_augmented);
    239307    free(solver->rScale);
    240308    free(solver->cScale);
    241309    free(solver->rwork);
    242310    free(solver->zero_vector);
     311    free(solver->iwork_array);
     312    free(solver->dgesdd_work);
     313    free(solver->dgesdd_iwork);
    243314    free(solver);
    244315    block->solver = 0;
    245316}