numerics/src/GenericMechanical/GMPReduced.c
/*
 * The equalities are assamblate in one block.
 *
 *0=(Me_1 Me_2)(Re Ri)' + Qe
 *Vi=(Mi_1 Mi_2)(Re Ri)' + Qi
 *
 *and GS.
 */
The equalities are assamblate in one block.

and GS.
The equalities are assamblate in one block . and GS .
void gmp_reduced_equality_solve(GenericMechanicalProblem* pInProblem, double *reaction, double *velocity, int * info, SolverOptions* options)
{
  SparseBlockStructuredMatrix* m = pInProblem->M->matrix1;
  int nbRow = m->blocksize0[m->blocknumber0 - 1];
  int nbCol = m->blocksize1[m->blocknumber1 - 1];
  int Me_size;
  int Mi_size;
  double * reducedProb = (double *)malloc(nbRow * nbCol * sizeof(double));
  double * Qreduced = (double *)malloc(nbRow * sizeof(double));
  double *Rreduced = (double *) calloc(nbCol, sizeof(double));
  double *Vreduced = (double *) calloc(nbRow, sizeof(double));
  _GMPReducedEquality(pInProblem, reducedProb, Qreduced, &Me_size, &Mi_size);
  if(Me_size == 0)
  {
    gmp_gauss_seidel(pInProblem, reaction, velocity, info, options);
    free(reducedProb);
    free(Qreduced);
    free(Rreduced);
    free(Vreduced);
    return;
  }
  listNumericsProblem * curProblem = 0;
  GenericMechanicalProblem * _pnumerics_GMP = genericMechanicalProblem_new();
  curProblem =  pInProblem->firstListElem;
  if(Me_size)
    gmp_add(_pnumerics_GMP, SICONOS_NUMERICS_PROBLEM_EQUALITY, Me_size);
  unsigned int curPos = 0;
  unsigned int curPosEq = 0;
  unsigned int curPosInq = Me_size;
  while(curProblem)
  {
    unsigned int curSize = curProblem->size;
    switch(curProblem->type)
    {
    case SICONOS_NUMERICS_PROBLEM_EQUALITY:
    {
      memcpy(Vreduced + curPosEq, velocity + curPos, curSize * sizeof(double));
      memcpy(Rreduced + curPosEq, reaction + curPos, curSize * sizeof(double));
      curPosEq += curSize;
      curPos += curSize;
      break;
    }
    case SICONOS_NUMERICS_PROBLEM_LCP:
    {
      memcpy(Vreduced + curPosInq, velocity + curPos, curSize * sizeof(double));
      memcpy(Rreduced + curPosInq, reaction + curPos, curSize * sizeof(double));
      curPosInq += curSize;
      curPos += curSize;
      gmp_add(_pnumerics_GMP, curProblem->type, curProblem->size);
      break;
    }
    case SICONOS_NUMERICS_PROBLEM_FC3D:
    {
      memcpy(Vreduced + curPosInq, velocity + curPos, curSize * sizeof(double));
      memcpy(Rreduced + curPosInq, reaction + curPos, curSize * sizeof(double));
      curPosInq += curSize;
      curPos += curSize;
      FrictionContactProblem* pFC3D = (FrictionContactProblem*)gmp_add(_pnumerics_GMP, curProblem->type, curProblem->size);
      *(pFC3D->mu) = *(((FrictionContactProblem*)curProblem->problem)->mu);
      break;
    }
    default:
      printf("GMPReduced  buildReducedGMP: problemType unknown: %d . \n", curProblem->type);
    }
    curProblem = curProblem->nextProblem;
  }
#ifdef GMP_DEBUG_GMPREDUCED_SOLVE
#endif
  NumericsMatrix numM;
  numM.storageType = 0;
  numM.matrix0 = reducedProb;
  numM.matrix1 = 0;
  numM.size0 = nbRow;
  numM.size1 = nbCol;
  _pnumerics_GMP->M = &numM;
  _pnumerics_GMP->q = Qreduced;
  gmp_gauss_seidel(_pnumerics_GMP, Rreduced, Vreduced, info, options);
#ifdef GMP_DEBUG_GMPREDUCED_SOLVE
  if(*info)
  {
    printf("\nGMPREduced2 failed!\n");
  }
  else
  {
    printf("\nGMPREduced2 succed!\n");
  }
#endif
  if(!*info)
  {
    gmp_reduced_convert_solution(pInProblem, reaction, velocity, Rreduced, Rreduced + Me_size, Vreduced + Me_size);
#ifdef GMP_DEBUG_GMPREDUCED_SOLVE
#endif
    double err;
    int tolViolate = gmp_compute_error(pInProblem, reaction, velocity, options->dparam[SICONOS_DPARAM_TOL], options, &err);
    if(tolViolate)
    {
      printf("GMPReduced2, warnning, reduced problem solved, but error of initial probleme violated tol = %e, err= %e\n", options->dparam[SICONOS_DPARAM_TOL], err);
    }
  }
#ifdef GMP_DEBUG_GMPREDUCED_SOLVE
#endif
  free(Rreduced);
  free(Vreduced);
  genericMechanicalProblem_free(_pnumerics_GMP, NUMERICS_GMP_FREE_GMP);
  free(Qreduced);
  free(reducedProb);
}
void gmp_reduced_equality_solve ( GenericMechanicalProblem * pInProblem , double * reaction , double * velocity , int * info , SolverOptions * options ) { SparseBlockStructuredMatrix * m = pInProblem -> M -> matrix1 ; int nbRow = m -> blocksize0 [ m -> blocknumber0 - 1 ] ; int nbCol = m -> blocksize1 [ m -> blocknumber1 - 1 ] ; int Me_size ; int Mi_size ; double * reducedProb = ( double * ) malloc ( nbRow * nbCol * sizeof ( double ) ) ; double * Qreduced = ( double * ) malloc ( nbRow * sizeof ( double ) ) ; double * Rreduced = ( double * ) calloc ( nbCol , sizeof ( double ) ) ; double * Vreduced = ( double * ) calloc ( nbRow , sizeof ( double ) ) ; _GMPReducedEquality ( pInProblem , reducedProb , Qreduced , & Me_size , & Mi_size ) ; if ( Me_size == 0 ) { gmp_gauss_seidel ( pInProblem , reaction , velocity , info , options ) ; free ( reducedProb ) ; free ( Qreduced ) ; free ( Rreduced ) ; free ( Vreduced ) ; return ; } listNumericsProblem * curProblem = 0 ; GenericMechanicalProblem * _pnumerics_GMP = genericMechanicalProblem_new ( ) ; curProblem = pInProblem -> firstListElem ; if ( Me_size ) gmp_add ( _pnumerics_GMP , SICONOS_NUMERICS_PROBLEM_EQUALITY , Me_size ) ; unsigned int curPos = 0 ; unsigned int curPosEq = 0 ; unsigned int curPosInq = Me_size ; while ( curProblem ) { unsigned int curSize = curProblem -> size ; switch ( curProblem -> type ) { case SICONOS_NUMERICS_PROBLEM_EQUALITY : { memcpy ( Vreduced + curPosEq , velocity + curPos , curSize * sizeof ( double ) ) ; memcpy ( Rreduced + curPosEq , reaction + curPos , curSize * sizeof ( double ) ) ; curPosEq += curSize ; curPos += curSize ; break ; } case SICONOS_NUMERICS_PROBLEM_LCP : { memcpy ( Vreduced + curPosInq , velocity + curPos , curSize * sizeof ( double ) ) ; memcpy ( Rreduced + curPosInq , reaction + curPos , curSize * sizeof ( double ) ) ; curPosInq += curSize ; curPos += curSize ; gmp_add ( _pnumerics_GMP , curProblem -> type , curProblem -> size ) ; break ; } case SICONOS_NUMERICS_PROBLEM_FC3D : { memcpy ( Vreduced + curPosInq , velocity + curPos , curSize * sizeof ( double ) ) ; memcpy ( Rreduced + curPosInq , reaction + curPos , curSize * sizeof ( double ) ) ; curPosInq += curSize ; curPos += curSize ; FrictionContactProblem * pFC3D = ( FrictionContactProblem * ) gmp_add ( _pnumerics_GMP , curProblem -> type , curProblem -> size ) ; * ( pFC3D -> mu ) = * ( ( ( FrictionContactProblem * ) curProblem -> problem ) -> mu ) ; break ; } default : printf ( " \n " , curProblem -> type ) ; } curProblem = curProblem -> nextProblem ; } #ifdef GMP_DEBUG_GMPREDUCED_SOLVE #endif NumericsMatrix numM ; numM . storageType = 0 ; numM . matrix0 = reducedProb ; numM . matrix1 = 0 ; numM . size0 = nbRow ; numM . size1 = nbCol ; _pnumerics_GMP -> M = & numM ; _pnumerics_GMP -> q = Qreduced ; gmp_gauss_seidel ( _pnumerics_GMP , Rreduced , Vreduced , info , options ) ; #ifdef GMP_DEBUG_GMPREDUCED_SOLVE if ( * info ) { printf ( " \n \n " ) ; } else { printf ( " \n \n " ) ; } #endif if ( ! * info ) { gmp_reduced_convert_solution ( pInProblem , reaction , velocity , Rreduced , Rreduced + Me_size , Vreduced + Me_size ) ; #ifdef GMP_DEBUG_GMPREDUCED_SOLVE #endif double err ; int tolViolate = gmp_compute_error ( pInProblem , reaction , velocity , options -> dparam [ SICONOS_DPARAM_TOL ] , options , & err ) ; if ( tolViolate ) { printf ( " \n " , options -> dparam [ SICONOS_DPARAM_TOL ] , err ) ; } } #ifdef GMP_DEBUG_GMPREDUCED_SOLVE #endif free ( Rreduced ) ; free ( Vreduced ) ; genericMechanicalProblem_free ( _pnumerics_GMP , NUMERICS_GMP_FREE_GMP ) ; free ( Qreduced ) ; free ( reducedProb ) ; }
The equalities are assamblate in one block.
The equalities are assamblate in one block .
//  printDenseMatrice("newPrb",file,reducedProb,nbRow,nbCol); //  printDenseMatrice("newQ",file,Qreduced,nbRow,1); //    printDenseMatrice("R",file,Rreduced,nbRow,1); //    printDenseMatrice("V",file,Vreduced,nbRow,1); //    printDenseMatrice("R2",file,reaction,nbRow,1); //    printDenseMatrice("V2",file,velocity,nbRow,1); //    fclose(file);