Changeset 9641


Ignore:
Timestamp:
Feb 16, 2017 9:29:18 AM (3 years ago)
Author:
Christian Andersson
Message:

Minor improvements to the sparse linear solver. Related to ticket:4815

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

Legend:

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

    r9597 r9641  
    2424#include "jmi_log.h"
    2525#include <stdint.h>
     26
    2627
    2728#define SMALL 1e-15
     
    175176       
    176177        for (i = 0; i < n_x-1; i++) {
    177             set = solver->dependent_set[(n_x-1)*n_x+i];
     178            set = (int)solver->dependent_set[(n_x-1)*n_x+i];
    178179            nbr_active = 0;
    179180           
     
    620621/* C (dense) = -A (sparse)*B (sparse) */
    621622int jmi_linear_solver_sparse_multiply(const jmi_matrix_sparse_csc_t *A, const jmi_matrix_sparse_csc_t *B, double *C) {
    622     jmi_int_t B_col, B_p, p;
     623    jmi_int_t B_col;
    623624
    624625    for (B_col = 0; B_col < B->nbr_cols; B_col++) {
    625         jmi_int_t col_ind = A->nbr_rows*B_col;
    626        
    627         for (B_p = B->col_ptrs[B_col]; B_p < B->col_ptrs[B_col+1]; B_p++) {
    628             double val = B->x[B_p];
    629             jmi_int_t col = B->row_ind[B_p];
    630            
    631             for (p = A->col_ptrs[col]; p < A->col_ptrs[col+1]; p++) {
    632                 C[col_ind + A->row_ind[p]] -=  A->x[p]*val;
    633             }
     626        jmi_linear_solver_sparse_multiply_column(A, B, B_col, C);
     627    }
     628   
     629    return 0;
     630}
     631
     632/* C (dense) = -A (sparse)*B(:,col) (sparse) */
     633int jmi_linear_solver_sparse_multiply_column(const jmi_matrix_sparse_csc_t *A, const jmi_matrix_sparse_csc_t *B, jmi_int_t B_col, double *C) {
     634    jmi_int_t B_p, p;
     635    jmi_int_t col_ind = A->nbr_rows*B_col;
     636
     637    for (B_p = B->col_ptrs[B_col]; B_p < B->col_ptrs[B_col + 1]; B_p++) {
     638        double val = B->x[B_p];
     639        jmi_int_t col = B->row_ind[B_p];
     640
     641        for (p = A->col_ptrs[col]; p < A->col_ptrs[col + 1]; p++) {
     642            C[col_ind + A->row_ind[p]] -= A->x[p] * val;
    634643        }
    635644    }
     
    665674    /* Compute L^(-1) A12 */
    666675    if (Jsp->L != NULL ) {
    667         jmi_int_t col, i, k = 0;
    668        
    669         for (col = 0; col < Jsp->A12->nbr_cols; col++) {
    670 
    671             jmi_linear_solver_sparse_backsolve(Jsp->L, Jsp->A12, Jsp->nz_patterns[col], Jsp->nz_sizes[col], col, Jsp->work_x);
     676        jmi_int_t col;
     677
     678        memset(block->J->data, 0, sizeof(double)*block->n*block->n);
     679       
     680        {
     681            jmi_int_t tid = 0;
     682            double *work;
     683            work = Jsp->work_x[tid];
     684
     685            for (col = 0; col < Jsp->A12->nbr_cols; col++) {
     686                jmi_int_t i;
     687                jmi_int_t offset = Jsp->nz_offsets[col];
     688
     689                jmi_linear_solver_sparse_backsolve(Jsp->L, Jsp->A12, Jsp->nz_patterns[col], Jsp->nz_sizes[col], col, work);
    672690           
    673             for (i = 0; i < Jsp->nz_sizes[col]; i++) {
    674                 Jsp->M1->x[k] = Jsp->work_x[Jsp->nz_patterns[col][i]];
    675                 Jsp->work_x[Jsp->nz_patterns[col][i]] = 0.0; /* Reset work vector */
    676                 k = k + 1;
    677             }
    678         }
    679         memset(block->J->data, 0, sizeof(double)*block->n*block->n);
    680        
    681         /* Compute A21L^(-1)A12 */
    682         jmi_linear_solver_sparse_multiply(Jsp->A21, Jsp->M1, block->J->data);
     691                for (i = 0; i < Jsp->nz_sizes[col]; i++) {
     692                    Jsp->M1->x[offset+i] = work[Jsp->nz_patterns[col][i]];
     693                    work[Jsp->nz_patterns[col][i]] = 0.0; /* Reset work vector */
     694                }
     695
     696                /* Compute A21L^(-1)A12 */
     697                jmi_linear_solver_sparse_multiply_column(Jsp->A21, Jsp->M1, col, block->J->data);
     698            }
     699        }
     700       
     701       
    683702        /* Compute A22 - A21L^(-1)A12 */
    684703        jmi_linear_solver_sparse_add_inplace(Jsp->A22, block->J->data);
     
    722741        jmi_int_t *work;
    723742        jmi_int_t max_dim = Jsp->L->nbr_cols > Jsp->A22->nbr_cols ? Jsp->L->nbr_cols : Jsp->A22->nbr_cols;
     743        jmi_int_t max_threads = 1;
    724744       
    725745        work_nz_pattern   = (jmi_int_t*)calloc(Jsp->L->nbr_cols+1, sizeof(jmi_int_t));
    726746        work              = (jmi_int_t*)calloc(Jsp->L->nbr_cols, sizeof(jmi_int_t));
    727747       
     748        Jsp->nz_offsets   = (jmi_int_t*)calloc(Jsp->A12->nbr_cols, sizeof(jmi_int_t));
    728749        Jsp->nz_sizes     = (jmi_int_t*)calloc(Jsp->A12->nbr_cols, sizeof(jmi_int_t));
    729750        Jsp->nz_patterns  = (jmi_int_t**)calloc(Jsp->A12->nbr_cols, sizeof(jmi_int_t*));
    730         Jsp->work_x       = (double*)calloc(max_dim, sizeof(double));
     751
     752       
     753        /* Allocate work arrays for the different threads */
     754        Jsp->max_threads = max_threads;
     755        Jsp->work_x       = (double**)calloc(max_threads, sizeof(double*));
     756        for (col = 0; col < max_threads; col++) {
     757            Jsp->work_x[col] = (double*)calloc(max_dim, sizeof(double));
     758        }
     759
     760        Jsp->nz_offsets[col] = 0;
    731761   
    732762        /* Compute the sparsity structure of L^(-1) A12 */
     
    739769            Jsp->nz_patterns[col] = (jmi_int_t*)calloc(Jsp->nz_sizes[col], sizeof(jmi_int_t));
    740770            for (i = 0; i <  Jsp->nz_sizes[col]; i++) { Jsp->nz_patterns[col][i] = work_nz_pattern[i]; }
     771
     772            if (col < Jsp->A12->nbr_cols - 1) {
     773                Jsp->nz_offsets[col+1] = Jsp->nz_offsets[col] + Jsp->nz_sizes[col];
     774            }
    741775           
    742776            nzmax += Jsp->nz_sizes[col];
     
    763797        return -1;
    764798    }
     799   
     800    if(block->callbacks->log_options.log_level >= 4) {
     801        jmi_log_node_t node;
     802        node = jmi_log_enter_fmt(block->log, logInfo, "LinearSparsity", "Sparsity information in <block:%s>", block->label);
     803       
     804        if (Jsp->L != NULL)
     805            jmi_log_fmt(block->log, node, logInfo, "Torn matrix L <numberOfColumns: %d> <numberOfRows: %d> <nonZeroElements: %d>", Jsp->L->nbr_cols, Jsp->L->nbr_rows, Jsp->L->nnz);
     806        if (Jsp->A12 != NULL)
     807            jmi_log_fmt(block->log, node, logInfo, "Torn matrix A12 <numberOfColumns: %d> <numberOfRows: %d> <nonZeroElements: %d>", Jsp->A12->nbr_cols, Jsp->A12->nbr_rows, Jsp->A12->nnz);
     808        if (Jsp->A21 != NULL)
     809            jmi_log_fmt(block->log, node, logInfo, "Torn matrix A21 <numberOfColumns: %d> <numberOfRows: %d> <nonZeroElements: %d>", Jsp->A21->nbr_cols, Jsp->A21->nbr_rows, Jsp->A21->nnz);
     810        if (Jsp->A22 != NULL)
     811            jmi_log_fmt(block->log, node, logInfo, "Torn matrix A22 <numberOfColumns: %d> <numberOfRows: %d> <nonZeroElements: %d>", Jsp->A22->nbr_cols, Jsp->A22->nbr_rows, Jsp->A22->nnz);
     812        if (Jsp->M1 != NULL)
     813            jmi_log_fmt(block->log, node, logInfo, "Torn matrix L^(-1)A12 <numberOfColumns: %d> <numberOfRows: %d> <nonZeroElements: %d>", Jsp->M1->nbr_cols, Jsp->M1->nbr_rows, Jsp->M1->nnz);
     814       
     815        jmi_log_leave(block->log, node);
     816    }
    765817
    766818    return 0;
     
    780832        }
    781833        if (Jsp->nz_sizes != NULL)  { free(Jsp->nz_sizes); }
    782         if (Jsp->work_x != NULL)    { free(Jsp->work_x); }
     834        if (Jsp->work_x != NULL)    {
     835            int i;
     836            for (i = 0; i < Jsp->max_threads; i++) {
     837                if (Jsp->work_x[i] != NULL) { free(Jsp->work_x[i]); }
     838            }
     839            free(Jsp->work_x);
     840        }
    783841
    784842        if (Jsp->L != NULL)   { jmi_linear_solver_delete_sparse_matrix(Jsp->L);   }
     
    889947       
    890948        if((block->callbacks->log_options.log_level >= 6)) {
    891             jmi_log_node_t node;
     949            jmi_log_node_t node;
    892950            node = jmi_log_enter_fmt(block->log, logInfo, "LinearSaveState", "Saving the Linear state in <block:%s>", block->label);
    893951            jmi_log_reals(block->log, node, logInfo, "ivs", block->last_accepted_x, block->n);
  • trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.h

    r9569 r9641  
    5151int jmi_linear_solver_solve(jmi_block_solver_t* block);
    5252
    53 int jmi_linear_solver_evaluate_jacobian(jmi_block_solver_t* block, jmi_real_t* jacobian);
    54 
    5553int jmi_linear_solver_evaluate_jacobian_factorization(jmi_block_solver_t* block, jmi_real_t* factorization);
    5654
     
    7068int jmi_linear_solver_sparse_multiply(const jmi_matrix_sparse_csc_t *A, const jmi_matrix_sparse_csc_t *B, double *C);
    7169
     70/** \brief Computes C(:,col) (dense) = -A (sparse)*B(:,col) (sparse) */
     71int jmi_linear_solver_sparse_multiply_column(const jmi_matrix_sparse_csc_t *A, const jmi_matrix_sparse_csc_t *B, jmi_int_t B_col, double *C);
     72
    7273/** \brief Computes C (dense) += A (sparse) */
    7374int jmi_linear_solver_sparse_add_inplace(const jmi_matrix_sparse_csc_t *A, double *C);
     
    8485    jmi_int_t** nz_patterns;
    8586    jmi_int_t* nz_sizes;
    86     double *work_x;
     87    jmi_int_t* nz_offsets;
     88    double **work_x;
     89    jmi_int_t max_threads;
    8790};
    8891
Note: See TracChangeset for help on using the changeset viewer.