#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: mlCheck.c,v 1.2 2000/01/10 03:20:33 knepley Exp $";
#endif
/*
   Defines error checking routines for the multilevel preconditioner
*/
#include "src/sles/pc/pcimpl.h"    /*I "pc.h" I*/
#include "ml.h"

#undef  __FUNCT__
#define __FUNCT__ "PCValidQ_Multilevel"
/*
  PCValidQ_Multilevel - Validates the ML data structure

  Collective on PC

  Input Parameter:
. pc - The PC

  Level: advanced

.keywords PC, Valid
.seealso PCDebug_Multilevel()
*/
int PCValidQ_Multilevel(PC pc)
{
  PC_Multilevel *ml = (PC_Multilevel *) pc->data;
  int            numRows, numCols;
  int            level, part;
  int            ierr;

  PetscFunctionBegin;
  if (pc->setupcalled < 2)
    PetscFunctionReturn(1);

  PetscValidPointer(ml);
  /* Check dimensions */
  if (ml->locB != PETSC_NULL) {
    ierr = MatGetSize(ml->locB, &numRows, &numCols);                                                     CHKERRQ(ierr);
    if ((ml->numRows != numRows) || (ml->numCols != numCols)) {
      PetscLogInfo(pc, "PCValidQ_Multilevel: Invalid dimensions (%d,%d) for factorization", numRows, numCols);
      PetscFunctionReturn(1);
    }
  }
  if (ml->numLevels < 0) {
    PetscLogInfo(pc, "PCValidQ_Multilevel: Invalid number of levels %d for factorization", ml->numLevels);
    PetscFunctionReturn(1);
  }
  /* Check thresholds */
  if (ml->QRthresh < 0) {
    PetscLogInfo(pc, "PCValidQ_Multilevel: Invalid threshold %d for final QR factorization", ml->QRthresh);
    PetscFunctionReturn(1);
  }
  if (ml->zeroTol < PETSC_MACHINE_EPSILON) {
    PetscLogInfo(pc, "PCValidQ_Multilevel: Numeric tolerance %g less than machine epsilon", ml->zeroTol);
    PetscFunctionReturn(1);
  }

  /* Check meshes */
  if (ml->numMeshes != ml->numLevels + 1)
    PetscLogInfo(pc, "PCValidQ_Multilevel: Invalid number %d of hierarchical meshes", ml->numMeshes);
  PetscValidPointer(ml->numElements);
  PetscValidPointer(ml->numVertices);
  PetscValidPointer(ml->meshes);
  for(level = 0; level <= ml->numLevels; level++) {
    PetscValidPointer(ml->meshes[level]);
    PetscValidPointer(ml->meshes[level][MESH_OFFSETS]);
    PetscValidPointer(ml->meshes[level][MESH_ADJ]);
    PetscValidPointer(ml->meshes[level][MESH_ELEM]);
  }

  /* Bail out if no factorization was done */
  if (ml->numLevels == 0)
    PetscFunctionReturn(0);

  /* Check numbering */
  PetscValidPointer(ml->numPartitions);
  PetscValidPointer(ml->numPartitionCols);
  PetscValidPointer(ml->colPartition);
  PetscValidPointer(ml->numPartitionRows);
  PetscValidPointer(ml->rowPartition);
  for(level = 0; level < ml->numLevels; level++) {
    if (ml->numPartitions[level] == 0)
      continue;
    PetscValidPointer(ml->numPartitionCols[level]);
    PetscValidPointer(ml->colPartition[level]);
    PetscValidPointer(ml->numPartitionRows[level]);
    PetscValidPointer(ml->rowPartition[level]);
    PetscValidPointer(ml->rowPartition[level][PART_ROW_INT]);
    PetscValidPointer(ml->rowPartition[level][PART_ROW_BD]);
    PetscValidPointer(ml->rowPartition[level][PART_ROW_RES]);
    if (ml->numPartitionRows[level][ml->numPartitions[level]] > 0)
      {PetscValidPointer(ml->rowPartition[level][PART_ROW_BD][0]);}
    if (ml->numPartitionRows[level][ml->numPartitions[level]+1] > 0)
      {PetscValidPointer(ml->rowPartition[level][PART_ROW_RES][0]);}
    for(part = 0; part < ml->numPartitions[level]; part++) {
      if (ml->numPartitionCols[level][part] > 0)
        {PetscValidPointer(ml->colPartition[level][part]);}
      if (ml->numPartitionRows[level][part] > 0)
        {PetscValidPointer(ml->rowPartition[level][PART_ROW_INT][part]);}
    }
  }

  /* Check factorization */
  PetscValidPointer(ml->factors);
  for(level = 0; level < ml->numLevels; level++) {
    if (ml->numPartitions[level] == 0)
      continue;
    PetscValidPointer(ml->factors[level]);
    for(part = 0; part < ml->numPartitions[level]; part++) {
      numRows = ml->numPartitionRows[level][part];
      numCols = ml->numPartitionCols[level][part];
      PetscValidPointer(ml->factors[level][part]);
      if (numRows > 0) {
        PetscValidPointer(ml->factors[level][part][FACT_U]);
        if (PCMultiLevelDoQR_Private(pc, numRows, numCols) == PETSC_TRUE) {
          PetscValidPointer(ml->factors[level][part][FACT_QR]);
          PetscValidPointer(ml->factors[level][part][FACT_TAU]);
        }
      }
      if (numCols > 0) {
        PetscValidPointer(ml->factors[level][part][FACT_DINV]);
        PetscValidPointer(ml->factors[level][part][FACT_V]);
      }
    }
  }

  /* Check boundary gradients */
  PetscValidPointer(ml->grads);
  PetscValidPointer(ml->bdReduceVecs);
  PetscValidPointer(ml->colReduceVecs);
  PetscValidPointer(ml->colReduceVecs2);
  PetscValidPointer(ml->interiorWork);
  PetscValidPointer(ml->interiorWork2);
  for(level = 0; level < ml->numLevels; level++) {
    if (ml->numPartitions[level] == 0)
      continue;
    PetscValidHeaderSpecific(ml->grads[level],          MAT_COOKIE);
    PetscValidHeaderSpecific(ml->bdReduceVecs[level],   VEC_COOKIE);
    PetscValidHeaderSpecific(ml->colReduceVecs[level],  VEC_COOKIE);
    PetscValidHeaderSpecific(ml->colReduceVecs2[level], VEC_COOKIE);
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCDebugPrintMat_Multilevel_Private"
int PCDebugPrintMat_Multilevel_Private(PC pc, VarOrdering rowOrder, VarOrdering colOrder, int (*applyM)(PC, GVec, GVec))
{
  PC_Multilevel *ml   = (PC_Multilevel *) pc->data;
  Grid           grid = ml->grid;
  PetscScalar    zero = 0.0;
  GVec           rowVec,    colVec;
  PetscScalar   *rowArray, *colArray;
  int            rank;
  int            tempLevels, colVars, locRowVars;
  int           *firstCol;
  int            row, col;
  PetscTruth     opt;
  int            ierr;

  PetscFunctionBegin;
  /* Setup storage */
  ierr = MPI_Comm_rank(pc->comm, &rank);                                                                 CHKERRQ(ierr);
  ierr = GVecCreateRectangular(grid, rowOrder, &rowVec);                                                 CHKERRQ(ierr);
  ierr = GVecCreateRectangular(grid, colOrder, &colVec);                                                 CHKERRQ(ierr);
  ierr = VecGetArray(rowVec, &rowArray);                                                                 CHKERRQ(ierr);
  ierr = VecGetArray(colVec, &colArray);                                                                 CHKERRQ(ierr);

  /* Allow only the first k levels to be used */
  tempLevels = ml->numLevels;
  ierr = PetscOptionsGetInt(pc->prefix, "-pc_ml_debug_level", &tempLevels, &opt);                         CHKERRQ(ierr);

  /* Print out the matrix M */
  ierr = VarOrderingGetSize(colOrder, &colVars);                                                         CHKERRQ(ierr);
  ierr = VarOrderingGetLocalSize(rowOrder, &locRowVars);                                                 CHKERRQ(ierr);
  ierr = VarOrderingGetFirst(colOrder, &firstCol);                                                       CHKERRQ(ierr);
  for(col = 0; col < colVars; col++)
  {
    /* Start with e_col */
    ierr = VecSet(&zero, colVec);                                                                        CHKERRQ(ierr);
    if ((col >= firstCol[rank]) && (col < firstCol[rank+1]))
      colArray[col-firstCol[rank]] = 1.0;
    /* M_col = M e_col */
    ierr = (*applyM)(pc, colVec, rowVec);                                                                CHKERRQ(ierr);
    /* Print out as a row so we get the transpose */
    for(row = 0; row < locRowVars; row++)
      PetscSynchronizedPrintf(pc->comm, "%g ", rowArray[row]);
    PetscSynchronizedFlush(pc->comm);
    PetscPrintf(pc->comm, "\n");
  }

  /* Cleanup stroage */
  ierr = VecRestoreArray(rowVec, &rowArray);                                                             CHKERRQ(ierr);
  ierr = VecRestoreArray(colVec, &colArray);                                                             CHKERRQ(ierr);
  ierr = GVecDestroy(rowVec);                                                                            CHKERRQ(ierr);
  ierr = GVecDestroy(colVec);                                                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCDebugPrintMat_Multilevel"
/*
  PCDebugPrintMat_Multilevel - Prints out the dense form of the level matrices.

  Collective on PC

  Input Parameter:
. pc - The PC

  Level: advanced

.keywords PC, debug
.seealso PCDebug_Multilevel()
*/
int PCDebugPrintMat_Multilevel(PC pc, const char mat[])
{
  PC_Multilevel *ml = (PC_Multilevel *) pc->data;
  PetscTruth     isP, isB, isV, isD, isall;
  int            ierr;

  PetscFunctionBegin;
  ierr = PetscStrcasecmp(mat, "P",   &isP);                                                              CHKERRQ(ierr);
  ierr = PetscStrcasecmp(mat, "B",   &isB);                                                              CHKERRQ(ierr);
  ierr = PetscStrcasecmp(mat, "V",   &isV);                                                              CHKERRQ(ierr);
  ierr = PetscStrcasecmp(mat, "D",   &isD);                                                              CHKERRQ(ierr);
  ierr = PetscStrcasecmp(mat, "all", &isall);                                                            CHKERRQ(ierr);
  if (isP == PETSC_TRUE) {
    PetscPrintf(pc->comm, "P^T\n");
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->tOrder, ml->tOrder, PCMultiLevelApplyP);           CHKERRQ(ierr);
  } else if (isB == PETSC_TRUE) {
    PetscPrintf(pc->comm, "B^T\n");
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->tOrder, ml->sOrder, PCMultiLevelApplyGradient);    CHKERRQ(ierr);
  } else if (isV == PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD, "V^T\n");
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->sOrder, ml->sOrder, PCMultiLevelApplyV);           CHKERRQ(ierr);
  } else if (isD == PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD, "D^{-T}\n");
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->sOrder, ml->sOrder, PCMultiLevelApplyDInv);        CHKERRQ(ierr);
  } else if (isall == PETSC_TRUE) {
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->tOrder, ml->tOrder, PCMultiLevelApplyP);           CHKERRQ(ierr);
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->tOrder, ml->sOrder, PCMultiLevelApplyGradient);    CHKERRQ(ierr);
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->sOrder, ml->sOrder, PCMultiLevelApplyV);           CHKERRQ(ierr);
    ierr = PCDebugPrintMat_Multilevel_Private(pc, ml->sOrder, ml->sOrder, PCMultiLevelApplyDInv);        CHKERRQ(ierr);
  } else {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Unknown ML level matrix %s.", mat);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCDebug_Multilevel"
int PCDebug_Multilevel(PC pc)
{
  PC_Multilevel       *ml   = (PC_Multilevel *) pc->data;
  Grid                 grid = ml->grid;
  PetscScalar          zero = 0.0;
  GVec                 testVec,    testVec2,    testVec3;
  PetscScalar         *testArray, *testArray2, *testArray3;
  GVec                 singVec;
  PetscScalar         *singArray;
  PetscScalar          norm;
  char                 printMat[256];
  int                  rank;
  int                  range;
  int                  rowVars, colVars, locColVars;
  int                 *firstRow, *firstCol;
  int                  singCount,    singRow;
  int                  locSingCount, locSingRow;
  int                  row, col, locCol;
  PetscTruth           opt;
  int                  ierr;

  PetscFunctionBegin;
  /* Initialize testing */
  ierr = MPI_Comm_rank(pc->comm, &rank);                                                                 CHKERRQ(ierr);
  ierr = GVecCreateRectangular(grid, ml->sOrder, &testVec);                                              CHKERRQ(ierr);
  ierr = GVecDuplicate(testVec, &singVec);                                                               CHKERRQ(ierr);
  ierr = GVecCreateConstrained(grid, &testVec2);                                                         CHKERRQ(ierr);
  ierr = GVecCreateConstrained(grid, &testVec3);                                                         CHKERRQ(ierr);
  ierr = VecGetArray(testVec,  &testArray);                                                              CHKERRQ(ierr);
  ierr = VecGetArray(testVec2, &testArray2);                                                             CHKERRQ(ierr);
  ierr = VecGetArray(testVec3, &testArray3);                                                             CHKERRQ(ierr);
  ierr = VecGetArray(singVec,  &singArray);                                                              CHKERRQ(ierr);

  /* Print matrices */
  ierr = PetscOptionsGetString(pc->prefix, "-pc_ml_print_mat", printMat, 255, &opt);                      CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = PCDebugPrintMat_Multilevel(pc, printMat);                                                     CHKERRQ(ierr);
  }

  ierr = VarOrderingGetSize(ml->tOrder, &rowVars);                                                       CHKERRQ(ierr);
  ierr = VarOrderingGetSize(ml->sOrder, &colVars);                                                       CHKERRQ(ierr);
  ierr = VarOrderingGetLocalSize(ml->sOrder, &locColVars);                                               CHKERRQ(ierr);
  ierr = VarOrderingGetFirst(ml->tOrder, &firstRow);                                                     CHKERRQ(ierr);
  ierr = VarOrderingGetFirst(ml->sOrder, &firstCol);                                                     CHKERRQ(ierr);

  /* Check that B^T P = 0 for numNullSpace rows */
  for(row = 0, range = 0; row < rowVars; row++)
  {
    /* Start with e_row */
    ierr = VecSet(&zero, testVec2);                                                                      CHKERRQ(ierr);
    if ((row >= firstRow[rank]) && (row < firstRow[rank+1]))
      testArray2[row-firstRow[rank]] = 1.0;
    /* P e_row */
    ierr = PCMultiLevelApplyP(pc, testVec2, testVec2);                                                   CHKERRQ(ierr);
    /* B^T P e_row */
    ierr = PCMultiLevelApplyGradientTrans(pc, testVec2, testVec);                                        CHKERRQ(ierr);
    /* Check for nonzeros only in the range */
    ierr = VecNorm(testVec, NORM_2, &norm);                                                              CHKERRQ(ierr);
    if (norm > ml->zeroTol)
      range++;
  }
  if (range != ml->globalRank) {
    PetscLogInfo(pc, "PCDebug_Multilevel: P_2 is not a null space for B^T, %d != %d range vectors\n", range, ml->globalRank);
    PetscFunctionReturn(1);
  }

  /* Test to see that D^{-1} P^T_1 B Z = I */
  for(col = 0; col < colVars; col++)
  {
    /* Start with e_col */
    locCol = col - firstCol[rank];
    ierr = VecSet(&zero, testVec);                                                                       CHKERRQ(ierr);
    if ((col >= firstCol[rank]) && (col < firstCol[rank+1]))
      testArray[locCol] = 1.0;

    /* Z e_col */
    ierr = PCMultiLevelApplyV(pc, testVec, testVec);                                                     CHKERRQ(ierr);
    /* B Z e_col */
    ierr = PCMultiLevelApplyGradient(pc, testVec, testVec2);                                             CHKERRQ(ierr);
    /* P^T B Z e_col */
    ierr = PCMultiLevelApplyPTrans(pc, testVec2, testVec2);                                              CHKERRQ(ierr);
    /* Scatter to a column vector */
    ierr = VecScatterBegin(testVec2, testVec, INSERT_VALUES, SCATTER_FORWARD, ml->rangeScatter);         CHKERRQ(ierr);
    ierr = VecScatterEnd(testVec2, testVec, INSERT_VALUES, SCATTER_FORWARD, ml->rangeScatter);           CHKERRQ(ierr);
    /* D^{-1} P^T B Z e_col */
    ierr = PCMultiLevelApplyDInv(pc, testVec, testVec);                                                  CHKERRQ(ierr);
    /* Check the row */
    for(row = 0, locSingCount = 0, locSingRow = -1; row < locColVars; row++)
      if (testArray[row] > ml->zeroTol) {
        locSingCount++;
        locSingRow = row;
      }
    ierr = MPI_Allreduce(&locSingCount, &singCount, 1, MPI_INT, MPI_SUM, pc->comm);                      CHKERRQ(ierr);
    ierr = MPI_Allreduce(&locSingRow,   &singRow,   1, MPI_INT, MPI_MAX, pc->comm);                      CHKERRQ(ierr);
    if (singCount > 1) {
      PetscLogInfo(pc, "PCDebug_Multilevel: Invalid column %d in P^T B Z, %d nonzeros\n", col, singCount);
      PetscFunctionReturn(1);
    } else if (singCount == 0) {
      PetscLogInfo(pc, "PCDebug_Multilevel: B is rank deficient in column %d\n", col);
    }
    /* Check the singular value */
    if (locSingRow > -1)
    {
      singRow += firstCol[rank];
      if (singRow != col) {
        PetscLogInfo(pc, "PCDebug_Multilevel: Invalid ordering in P^T B Z, value in column %d and row %d\n", col, singRow);
        PetscFunctionReturn(1);
      }
      if (PetscAbsScalar(testArray[locSingRow] - 1.0) > ml->zeroTol) {
        PetscLogInfo(pc, "PCDebug_Multilevel: Invalid singular value in column %d\n", col);
        PetscFunctionReturn(1);
      }
    }
  }

  /* Cleanup testing */
  ierr = VecRestoreArray(testVec,  &testArray);                                                          CHKERRQ(ierr);
  ierr = VecRestoreArray(testVec2, &testArray2);                                                         CHKERRQ(ierr);
  ierr = VecRestoreArray(testVec3, &testArray3);                                                         CHKERRQ(ierr);
  ierr = VecRestoreArray(singVec,  &singArray);                                                          CHKERRQ(ierr);
  ierr = GVecDestroy(testVec);                                                                           CHKERRQ(ierr);
  ierr = GVecDestroy(testVec2);                                                                          CHKERRQ(ierr);
  ierr = GVecDestroy(testVec3);                                                                          CHKERRQ(ierr);
  ierr = GVecDestroy(singVec);                                                                           CHKERRQ(ierr);

  PetscFunctionReturn(0);
}