@@ -218,9 +218,9 @@ void Cvode::initialize()
218218 if (_idid < 0 )
219219 throw std::invalid_argument (" Cvode::initialize()" );
220220
221- // Use own jacobian matrix
222- _idid = CVDlsSetDenseJacFn (_cvodeMem, CV_JCallback);
223- if (_idid < 0 )
221+ // Use own jacobian matrix
222+ _idid = CVDlsSetDenseJacFn (_cvodeMem, CV_JCallback);
223+ if (_idid < 0 )
224224 throw std::invalid_argument (" CVode::initialize()" );
225225
226226 if (_dimZeroFunc)
@@ -570,7 +570,7 @@ int Cvode::calcJacobian(double t, long int N, N_Vector fHelp, N_Vector errorWeig
570570 try
571571 {
572572 int j,k,l;
573- double fnorm, minInc, *f_data, *fHelp_data , *errorWeight_data, h, srur, delta_inv;
573+ double fnorm, minInc, *f_data, *fHelp_data , *errorWeight_data, h, srur, delta_inv;
574574
575575 f_data = NV_DATA_S (fy);
576576 errorWeight_data = NV_DATA_S (errorWeight);
@@ -632,48 +632,48 @@ int Cvode::calcJacobian(double t, long int N, N_Vector fHelp, N_Vector errorWeig
632632 {
633633 j = _sparsePattern_leadindex[ii-1 ];
634634
635- }
636- while (j <_sparsePattern_leadindex[ii])
637- {
638- l = _sparsePattern_index[j];
639- k = l + ii * _dimSys;
640- // Jac->data[k] = (fHelp_data[l] - f_data[l])/_delta[l];
641- delta_inv = 1.0 /_delta[ii];
642- Jac->data [k] = (fHelp_data [l] - f_data[l])*delta_inv;
643- j++;
644- }
645- }
646- }
647- }
648-
649- /*
650- //Calculation of J without colouring
651- for (j = 0; j < N; j++)
652- {
635+ }
636+ while (j <_sparsePattern_leadindex[ii])
637+ {
638+ l = _sparsePattern_index[j];
639+ k = l + ii * _dimSys;
640+ // Jac->data[k] = (fHelp_data[l] - f_data[l])/_delta[l];
641+ delta_inv = 1.0 /_delta[ii];
642+ Jac->data [k] = (fHelp_data [l] - f_data[l])*delta_inv;
643+ j++;
644+ }
645+ }
646+ }
647+ }
648+
649+ /*
650+ //Calculation of J without colouring
651+ for (j = 0; j < N; j++)
652+ {
653+
653654
655+ //N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol);
654656
655- //N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol) ;
657+ _ysave[j] = y[j] ;
656658
657- _ysave [j] = y [j];
659+ y [j] += _delta [j];
658660
659- y[j] += _delta[j] ;
661+ calcFunction(t, y, fHelp_data) ;
660662
661- calcFunction(t, y, fHelp_data);
662-
663- y[j] = _ysave[j];
663+ y[j] = _ysave[j];
664664
665- delta_inv = 1.0/_delta[j];
666- N_VLinearSum(delta_inv, fHelp, -delta_inv, fy, jthCol);
665+ delta_inv = 1.0/_delta[j];
666+ N_VLinearSum(delta_inv, fHelp, -delta_inv, fy, jthCol);
667667
668- for(int i=0; i<_dimSys; ++i)
668+ for(int i=0; i<_dimSys; ++i)
669669 {
670670 Jac->data[i+j*_dimSys] = NV_Ith_S(jthCol,i);
671671 }
672672
673- //DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);
674- }
675- */
676-
673+ //DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);
674+ }
675+ */
676+
677677 } // workaround until exception can be catch from c- libraries
678678 catch (std::exception& ex)
679679 {
0 commit comments