numerics/src/tools/FischerBurmeister.c
/* Compute the jacobian of the mixed Fischer function */
Compute the jacobian of the mixed Fischer function
Compute the jacobian of the mixed Fischer function
void jacobianPhi_Mixed_FB(int sizeEq, int sizeIneq, double* restrict z, double* restrict F, double* restrict jacobianF, double* restrict jacobianPhiMatrix)
{
  assert(z != NULL);
  assert(F != NULL);
  assert(jacobianF != NULL);
  assert(jacobianPhiMatrix != NULL);
  cblas_dcopy((sizeEq + sizeIneq) * (sizeEq + sizeIneq), jacobianF, 1, jacobianPhiMatrix, 1);
  double ri, ai, bi;
  for(int i = sizeEq; i < sizeEq + sizeIneq; i++)
  {
    ri = sqrt(z[i] * z[i] +  F[i] * F[i]);
    if(ri > 0.0)
    {
      ai = z[i] / ri - 1.0;
      bi = F[i] / ri - 1.0;
    }
    else
    {
      ai = -1.0;
      bi = -1.0;
    }
    cblas_dscal(sizeEq + sizeIneq, bi, &jacobianPhiMatrix[i], sizeEq + sizeIneq);
    jacobianPhiMatrix[i * (sizeEq + sizeIneq) + i] += ai;
  }
}
void jacobianPhi_Mixed_FB ( int sizeEq , int sizeIneq , double * restrict z , double * restrict F , double * restrict jacobianF , double * restrict jacobianPhiMatrix ) { assert ( z != NULL ) ; assert ( F != NULL ) ; assert ( jacobianF != NULL ) ; assert ( jacobianPhiMatrix != NULL ) ; cblas_dcopy ( ( sizeEq + sizeIneq ) * ( sizeEq + sizeIneq ) , jacobianF , 1 , jacobianPhiMatrix , 1 ) ; double ri , ai , bi ; for ( int i = sizeEq ; i < sizeEq + sizeIneq ; i ++ ) { ri = sqrt ( z [ i ] * z [ i ] + F [ i ] * F [ i ] ) ; if ( ri > 0.0 ) { ai = z [ i ] / ri - 1.0 ; bi = F [ i ] / ri - 1.0 ; } else { ai = -1.0 ; bi = -1.0 ; } cblas_dscal ( sizeEq + sizeIneq , bi , & jacobianPhiMatrix [ i ] , sizeEq + sizeIneq ) ; jacobianPhiMatrix [ i * ( sizeEq + sizeIneq ) + i ] += ai ; } }
Compute the jacobian of the mixed Fischer function
Compute the jacobian of the mixed Fischer function
/* jacobianPhiMatrix is initialized with jacobianF */ /* jacobianPhiMatrix_ij = delta_ij*ai + bi * jacobianF_ij
       delta_ij being the Kronecker symbol
    */ /*        cblas_dscal(size, bi, &jacobianPhiMatrix[i], size);*/