Ticket #3879: linear_solution_check.patch

File linear_solution_check.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*/
     28static double jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block) {
     29    jmi_linear_solver_t* solver = block->solver;
     30    char norm = 'I';
     31    int N = block->n;
     32    double jnorm = 1.0;
     33    double rcond = 1.0;
     34    int info;
     35
     36    /* Compute infinity norm of J to be used with dgecon */
     37    jnorm = dlange_(&norm, &N, &N, solver->jacobian, &N, solver->rwork);
     38
     39    /* Compute reciprocal condition number */
     40    dgecon_(&norm, &N, solver->factorization, &N, &jnorm, &rcond, solver->rwork, solver->iwork_array,&info);
     41
     42    return 1.0/rcond;
     43}
     44
     45
    2646int jmi_linear_solver_new(jmi_linear_solver_t** solver_ptr, jmi_block_solver_t* block) {
    2747    jmi_linear_solver_t* solver= (jmi_linear_solver_t*)calloc(1,sizeof(jmi_linear_solver_t));
    2848/*    jmi_t* jmi = block->jmi;
     
    3858    solver->factorization = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    3959    solver->jacobian = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    4060    solver->jacobian_temp = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
     61    solver->jacobian_temp_augmented = (jmi_real_t*)calloc(n_x*(n_x+1),sizeof(jmi_real_t));
    4162    solver->singular_values = (jmi_real_t*)calloc(n_x,sizeof(jmi_real_t));
    4263    solver->rScale = (double*)calloc(n_x,sizeof(double));
    4364    solver->cScale = (double*)calloc(n_x,sizeof(double));
     
    4667   
    4768    solver->singular_jacobian = 0;
    4869    solver->iwork = 5*n_x;
     70    solver->iwork_array = (int*)calloc(solver->iwork,sizeof(int));
    4971    solver->rwork = (double*)calloc(solver->iwork,sizeof(double));
    5072    solver->zero_vector = (double*)calloc(n_x,sizeof(double));
     73   
     74    solver->dgesdd_lwork = 10*n_x;
     75    solver->dgesdd_work  = (double*)calloc(solver->dgesdd_lwork,sizeof(double));
     76    solver->dgesdd_iwork = (int*)calloc(8*n_x,sizeof(int));
    5177
    5278    for (i=0; i<n_x; i++) {
    5379        solver->zero_vector[i] = 0.0;
     
    164190    }
    165191   
    166192    if((block->callbacks->log_options.log_level >= 5)) {
     193        /* double rcond = jmi_calculate_jacobian_condition_number(block); */
    167194        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);     
     195        jmi_log_reals(block->log, destnode, logInfo, "b", block->res, block->n);
     196        /* jmi_log_reals(block->log, destnode, logInfo, "condition_number", &rcond, 1); */
    169197        jmi_log_leave(block->log, destnode);
    170198    }
    171199 
     
    174202    i = 1; /* One rhs to solve for */
    175203     
    176204    if (solver->singular_jacobian == 1){
     205        char jobz = 'N', ceps = 'P', csfmin = 'S';
     206        int j,rank_aug = 0, m_x = n_x + 1;
     207        double eps = dlamch_(&ceps);
     208        double sfmin = dlamch_(&csfmin);
     209        double threshold;
     210       
     211        memcpy(solver->jacobian_temp_augmented, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
     212        memcpy(&solver->jacobian_temp_augmented[n_x*n_x], block->res, n_x*sizeof(jmi_real_t));
    177213        /*
    178214         *   DGELSS - compute the minimum norm solution to      a real
    179215         *   linear least squares problem
     
    191227            return -1;
    192228        }
    193229       
     230        /* Verify solution */
     231        /*
     232         *       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
     233         *                          LWORK, IWORK, INFO )
     234         */
     235        dgesdd_(&jobz, &n_x, &m_x, solver->jacobian_temp_augmented, &n_x,
     236                solver->singular_values, NULL, &n_x, NULL, &n_x,
     237                solver->dgesdd_work, &solver->dgesdd_lwork, solver->dgesdd_iwork, &info);
     238               
     239        if(info != 0) {
     240            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);
     241            return -1;
     242        }
     243       
     244        threshold = (eps*solver->singular_values[0] > sfmin)? eps*solver->singular_values[0]: sfmin;
     245        for (j = 0; j < n_x; j++)
     246        {
     247            if (solver->singular_values[j] > threshold) {
     248                rank_aug += 1;
     249            } else {
     250                break;
     251            }
     252        }
     253       
     254        if (rank_aug != rank) {
     255            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);
     256            return -1;
     257        }
     258       
    194259        if(block->callbacks->log_options.log_level >= 5){
    195260            jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
    196261        }
     
    236301    free(solver->jacobian);
    237302    free(solver->singular_values);
    238303    free(solver->jacobian_temp);
     304    free(solver->jacobian_temp_augmented);
    239305    free(solver->rScale);
    240306    free(solver->cScale);
    241307    free(solver->rwork);
    242308    free(solver->zero_vector);
     309    free(solver->iwork_array);
     310    free(solver->dgesdd_work);
     311    free(solver->dgesdd_iwork);
    243312    free(solver);
    244313    block->solver = 0;
    245314}