matviewstruct.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2006-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: matviewstruct.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

int matviewstruct(MPI_Comm comm,Mat mat)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  PetscInt m,n;
  ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr);

  ierr = PetscViewerSetFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
  ierr = MatView(mat,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);

  if (m!=n) {
    PetscPrintf(PETSC_COMM_WORLD,"Matrix not square, cannot calculate bandwidth.\n");
    MagparFunctionLogReturn(0);
  }

/*
  ierr = MatView(mat,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
*/
/* very slow if running in parallel
  ierr = MatView(mat,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  MagparFunctionLogReturn(0);
*/
#if 0
  /* does not work 1 */
  Mat matloc;
  IS  isrow,iscol;
  ierr = MatGetSubMatrix(mat,isrow,iscol,n,MAT_INITIAL_MATRIX,&matloc);CHKERRQ(ierr);
#elif 0
  /* does not work 2 */
  PetscInt from,to;
  ierr = MatGetOwnershipRange(mat,&from,&to);CHKERRQ(ierr);

  IS  isrow,iscol;
  ierr = ISCreateStride(PETSC_COMM_WORLD,to-from,from,1,&isrow);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_WORLD,n,1,1,&iscol);CHKERRQ(ierr);

  Mat matloc;
  MatConvert(mat,MATSEQAIJ,MAT_INITIAL_MATRIX,&matloc);

/* this is the problem:
[1]PETSC ERROR: MatCreate_SeqAIJ() line 2864 in src/mat/impls/aij/seq/aij.c
[1]PETSC ERROR: Argument out of range!
[1]PETSC ERROR: Comm must be of size 1!
[1]PETSC ERROR: MatSetType() line 64 in src/mat/interface/matreg.c
[1]PETSC ERROR: MatConvert() line 3141 in src/mat/interface/matrix.c
*/

  ierr = ISDestroy(isrow);CHKERRQ(ierr);
  ierr = ISDestroy(iscol);CHKERRQ(ierr);
#elif 0
  /* does not work 3 */
  Mat matloc;
  ierr = MatGetLocalMat(mat,MAT_INITIAL_MATRIX,&matloc);CHKERRQ(ierr);
/*
  ierr = MatView(matloc,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  ierr = MatDestroy(matloc);CHKERRQ(ierr);
*/
#else
  PetscInt from,to;
  ierr = MatGetOwnershipRange(mat,&from,&to);CHKERRQ(ierr);

  IS  isrow,iscol;
  ierr = ISCreateStride(comm,to-from,from,1,&isrow);CHKERRQ(ierr);
  ierr = ISCreateStride(comm,n,1,1,&iscol);CHKERRQ(ierr);

  Mat *submatloc;
  ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&submatloc);CHKERRQ(ierr);

  ierr = ISDestroy(isrow);CHKERRQ(ierr);
  ierr = ISDestroy(iscol);CHKERRQ(ierr);

  Mat matloc;
  matloc=submatloc[0];
#endif

  /* calculate bandwidth */

  PetscInt *ia,*ja;
  PetscTruth done;
  ierr = MatGetRowIJ(matloc,0,PETSC_TRUE,PETSC_FALSE,&from,&ia,&ja,&done);CHKERRQ(ierr);

  /* assume that we get useful data back (not true on 2nd,...,nth processor)
  assert(done==PETSC_TRUE);
*/

  /* assume square matrix, no I-nodes used */
  int cmin,cmax,bw,bwmax;
  bwmax=0;
  for (int i=0;i<m;i++) {
    cmin=m;
    cmax=0;
    for (int j=0;j<(ia[i+1]-ia[i]);j++) {
      cmin=PetscMin(cmin,ja[ia[i]+j]);
      cmax=PetscMax(cmax,ja[ia[i]+j]);
    }
    bw=cmax-cmin+1;
    if (bw > bwmax) {
      bwmax=bw;
      n=i;
    }
  }
  PetscPrintf(comm,"max. bandwidth: %i in row %i (of %i)\n",bwmax,n,m);

  ierr = MatRestoreRowIJ(matloc,0,PETSC_TRUE,PETSC_FALSE,&from,&ia,&ja,&done);CHKERRQ(ierr);

#ifdef PETSC_HAVE_X11
  if (!rank) ierr = MatView(matloc,PETSC_VIEWER_DRAW_SELF);CHKERRQ(ierr);
#endif

  ierr = MatDestroyMatrices(1,&submatloc);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}