diff --git a/Makefile.am b/Makefile.am
index 0fc1bc545d9cd1794997881e790fbdf51b179ac0..1c48808ef6497e2f81adaaf992e4fcbd57c4ca17 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -60,6 +60,7 @@ libelpa@SUFFIX@_private_la_SOURCES = \
   src/elpa2/qr/elpa_pdgeqrf.F90 \
   src/elpa1/elpa1.F90 \
   src/elpa2/elpa2.F90 \
+  src/helpers/matrix_plot.F90 \
   src/elpa_c_interface.c \
   src/elpa_index.c
 
diff --git a/src/elpa1/elpa1_tridiag_template.X90 b/src/elpa1/elpa1_tridiag_template.X90
index 67a9717c8e6b6ca4fe9e6d9ff4cab10c5ddf396a..68b49d80e23f7a3e952a18c9700025d397a5d2c1 100644
--- a/src/elpa1/elpa1_tridiag_template.X90
+++ b/src/elpa1/elpa1_tridiag_template.X90
@@ -83,6 +83,15 @@
 !>
 #include "../general/sanity.X90"
 
+#undef SAVE_MATR
+#ifdef DOUBLE_PRECISION_REAL
+#define SAVE_MATR(name, iteration) \
+call prmat(na,useGpu,a_mat,a_dev,lda,matrixCols,nblk,my_prow,my_pcol,np_rows,np_cols,name,iteration) 
+#else
+#define SAVE_MATR(name, iteration)
+#endif
+
+
     subroutine tridiag_&
     &MATH_DATATYPE&
     &_&
@@ -92,6 +101,7 @@
       use iso_c_binding
       use precision
       use elpa_abstract_impl
+      use matrix_plot
       implicit none
       class(elpa_abstract_impl_t), intent(inout) :: obj
       integer(kind=ik), intent(in)                  :: na, lda, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols
@@ -446,6 +456,8 @@
           ! add tau after the end of actuall v_row, to be broadcasted with it
           v_row(l_rows+1) = tau(istep)
          endif !(my_pcol == pcol(istep, nblk, np_cols))
+         
+!          SAVE_MATR("HH vec stored", na - istep + 1)
 
 #ifdef WITH_MPI
          ! Broadcast the Householder Vector (and tau) along columns
@@ -1025,3 +1037,5 @@
     &MATH_DATATYPE&
     &_&
     &PRECISION
+
+       
\ No newline at end of file
diff --git a/src/helpers/matrix_plot.F90 b/src/helpers/matrix_plot.F90
new file mode 100644
index 0000000000000000000000000000000000000000..efc8c96f2e3380647d3dcdd1c887234628255ed3
--- /dev/null
+++ b/src/helpers/matrix_plot.F90
@@ -0,0 +1,80 @@
+#define  REAL_DATATYPE rk8
+
+! module for producing matrix traces, to be plotted by provided python ploter
+! currently the module is very simple and non-flexible
+! it is only ussable for printing the matrix A and possibly its counterpart A_DEV 
+! both are assumed to be in block-cyclic distribution
+! At the moment, the output works for double real only
+! To simplify thihngs, a conveninece macro (as follows) can be placed to .X90 template file:
+
+! #undef SAVE_MATR
+! #ifdef DOUBLE_PRECISION_REAL
+! #define SAVE_MATR(name, iteration) \
+! call prmat(na,useGpu,a_mat,a_dev,lda,matrixCols,nblk,my_prow,my_pcol,np_rows,np_cols,name,iteration) 
+! #else
+! #define SAVE_MATR(name, iteration)
+! #endif
+
+! traces are stored into directory "matrices", that has to be created
+
+module matrix_plot
+
+  contains
+  
+    subroutine prmat(na, useGpu, a_mat, a_dev, lda, matrixCols, nblk, my_prow, my_pcol, np_rows, np_cols, name, iteration)
+      use cuda_functions
+      use iso_c_binding
+      use precision
+      implicit none
+      integer, parameter :: out_unit=20
+      character(len = 1025) :: directory = "matrices"
+      character(len = 1024) :: filename
+      
+      character(len = *), intent(in)             :: name
+      integer(kind=ik), intent(in)                  :: na, lda, nblk, matrixCols, my_prow, my_pcol, np_rows, np_cols, iteration
+      real(kind=REAL_DATATYPE), intent(in)          :: a_mat(lda,matrixCols)
+      integer(kind=C_intptr_T), intent(in)          :: a_dev
+      logical, intent(in)                           :: useGPU
+
+      integer(kind=ik)                              :: row, col, mpi_rank
+      integer(kind=ik), save                        :: counter = 0
+      real(kind=REAL_DATATYPE)                      :: a_dev_helper(lda,matrixCols)
+      logical                                       :: successCUDA
+      integer(kind=c_size_t), parameter             :: size_of_datatype = size_of_double_real
+      
+      mpi_rank = np_rows * my_pcol + my_prow
+      
+      ! print a_mat
+      write(filename, "(A,A,I0.4,A,I0.2,A)"), trim(directory), "/a_mat-", counter, "-", mpi_rank, ".txt"    
+      write(*,*), trim(filename)
+      open(unit=out_unit, file=trim(filename), action="write",status="replace")
+      
+      write(out_unit, "(9I5)") na, nblk, lda, matrixCols, my_prow, my_pcol, np_rows, np_cols, iteration
+      write(out_unit, "(A)") name
+      do row = 1, lda
+          write(out_unit, *) a_mat(row, :)
+      end do
+      close(out_unit)    
+
+      ! print a_dev
+
+      if(useGpu) then
+        successCUDA = cuda_memcpy(loc(a_dev_helper(1,1)), a_dev, lda * matrixCols * size_of_datatype, cudaMemcpyDeviceToHost)
+
+        write(filename, "(A,A,I0.4,A,I0.2,A)"), trim(directory), "/a_dev-", counter, "-", mpi_rank, ".txt"    
+        write(*,*), trim(filename)
+        open(unit=out_unit, file=trim(filename), action="write",status="replace")
+        
+        write(out_unit, "(9I5)") na, nblk, lda, matrixCols, my_prow, my_pcol, np_rows, np_cols, iteration
+        write(out_unit, "(A)") name
+        do row = 1, lda
+            write(out_unit, *) a_dev_helper(row, :)
+        end do
+        close(out_unit)    
+      end if
+            
+      counter = counter + 1
+      
+    end subroutine
+
+end module
\ No newline at end of file
diff --git a/utils/matrix_plotter/defect_example/screenshot.png b/utils/matrix_plotter/defect_example/screenshot.png
new file mode 100644
index 0000000000000000000000000000000000000000..1526148ffc60bb1a3f5c43755c1ab61c263a0d04
Binary files /dev/null and b/utils/matrix_plotter/defect_example/screenshot.png differ
diff --git a/utils/matrix_plotter/defect_example/script.py b/utils/matrix_plotter/defect_example/script.py
new file mode 100644
index 0000000000000000000000000000000000000000..18a42d662f986b65f2aa54707b8580bca2ca7b15
--- /dev/null
+++ b/utils/matrix_plotter/defect_example/script.py
@@ -0,0 +1,46 @@
+# Sample script showing a use case of ELPA matrix output for debugging
+# In this case, we wanted to compare matrices for GPU and CPU variant of the code
+
+# This feature is only experimental and will probably change
+# The recomanded ussage is the following. Go to the directory of this script and run:
+# ipython -i ../plot.py 
+# than run within python the command:
+# run -i script.py
+# Now a window with the plot should be created. Rearange your windows such that you can 
+# see both console and plot and iterate using the iterators, change views or create new, etc.
+# If you want to modify this file and re-run it without closing the plotting window, comment
+# out the two following lines. If you do not want to reload the data, comment out creating of 
+# traces. 
+
+# data traces are too large for this example and are not included
+
+# create new matplotlib figure and pass it to the Plotter
+fig = plt.figure()
+pl = Plotter(fig)
+
+# load two traces, one from a run with GPU enabled, one without
+# note that this might take a lot of time, so if we have allready run this script and than only
+# changed something in the following, we might choose to comment out the following commands
+trace_ok = Trace('data_ok', "OK")
+trace_defect = Trace('data_defect', "defect")
+
+# descriptions are "metadata", they describe when and what has been outputed
+# we check that they do not differ for the two program runs
+assert(trace_ok._descriptions == trace_defect._descriptions)
+
+# we can then use one iterator to acces both data sets
+it = Iterator(trace_ok._descriptions)
+
+# we want to see matrix A both from the run we know was OK and that where was error
+a_mat_ok = Snapshot(trace_ok, "a_mat", it)
+a_mat_defect = Snapshot(trace_defect, "a_mat", it)
+
+# we also want to see where they differ, so that we can find out where the difference originates
+a_mat_diff = Difference(a_mat_ok, a_mat_defect)
+
+# we have created 3 views, lets group them into a list
+snapshot_set = [a_mat_ok, a_mat_defect, a_mat_diff]
+
+# and tell the Plotter to show them. Note, that all the plots were created to follow
+# the same iterator it, so we can change time point shown in all views by iterating it.
+pl.reset(snapshot_set)
diff --git a/utils/matrix_plotter/gpu_example/data.tgz b/utils/matrix_plotter/gpu_example/data.tgz
new file mode 100644
index 0000000000000000000000000000000000000000..a5586a38e0bfbaff74c0a9b539bbfafdab9bb598
Binary files /dev/null and b/utils/matrix_plotter/gpu_example/data.tgz differ
diff --git a/utils/matrix_plotter/gpu_example/screenshot.png b/utils/matrix_plotter/gpu_example/screenshot.png
new file mode 100644
index 0000000000000000000000000000000000000000..c678525fd08fc42348f24c98b19aa6e2effbbf0a
Binary files /dev/null and b/utils/matrix_plotter/gpu_example/screenshot.png differ
diff --git a/utils/matrix_plotter/gpu_example/script.py b/utils/matrix_plotter/gpu_example/script.py
new file mode 100644
index 0000000000000000000000000000000000000000..6611300551d613ca00b7e023eeb7dc766c2d7256
--- /dev/null
+++ b/utils/matrix_plotter/gpu_example/script.py
@@ -0,0 +1,69 @@
+# Sample script showing a use case of ELPA matrix output for debugging
+# In this case, we wanted to compare matrices for GPU and CPU variant of the code
+
+# This feature is only experimental and will probably change
+# The recomanded ussage is the following. Go to the directory of this script and run:
+# ipython -i ../plot.py 
+# than run within python the command:
+# run -i script.py
+# Now a window with the plot should be created. Rearange your windows such that you can 
+# see both console and plot and iterate using the iterators, change views or create new, etc.
+# If you want to modify this file and re-run it without closing the plotting window, comment
+# out the two following lines. If you do not want to reload the data, comment out creating of 
+# traces. 
+
+# To try this example untar the file data.tgz
+
+# create new matplotlib figure and pass it to the Plotter
+fig = plt.figure()
+pl = Plotter(fig)
+
+# load two traces, one from a run with GPU enabled, one without
+# note that this might take a lot of time, so if we have allready run this script and than only
+# changed something in the following, we might choose to comment out the following commands
+trace = Trace('data_gpu', "w gpu")
+trace_nongpu = Trace('data_nongpu', "nogpu")
+
+# descriptions are "metadata", they describe when and what has been outputed
+# we check that they do not differ for the two program runs
+assert(trace._descriptions == trace_nongpu._descriptions)
+
+# we can then use one iterator to acces both data sets
+it = Iterator(trace._descriptions)
+
+# we want to see two matrices from the run with GPU, we have both matrix on host and on device
+a_mat = Snapshot(trace, "a_mat", it)
+a_dev = Snapshot(trace, "a_dev", it)
+
+# we are also interested in difference between the two matrices
+a_mat_dev_diff = Difference(a_mat, a_dev)
+
+# from the second trace we can extract only matrix on host (no matrix on device)
+a_nongpu = Snapshot(trace_nongpu, "a_mat", it)
+
+# we have created 4 views, lets group them into a list
+snapshot_set = [a_mat, a_dev, a_mat_dev_diff, a_nongpu]
+
+# and tell the Plotter to show them. Note, that all the plots were created to follow
+# the same iterator it, so we can change time point shown in all views by iterating it.
+pl.reset(snapshot_set)
+
+# we can also create another iterator it_2 and the same shnapshot set, but binded to this new 
+# iterator it_2, to be able to switch between two time instances
+it_2 = Iterator(trace._descriptions)
+a_mat_2 = Snapshot(trace, "a_mat", it_2)
+a_dev_2 = Snapshot(trace, "a_dev", it_2)
+a_mat_dev_diff_2 = Difference(a_mat_2, a_dev_2)
+a_nongpu_2 = Snapshot(trace_nongpu, "a_mat", it_2)
+
+# we can then show this new set by calling pl.reset(snapshot_set_2)
+snapshot_set_2 = [a_mat_2, a_dev_2, a_mat_dev_diff_2, a_nongpu_2]
+
+# we have allready used Difference of two different matrices at the same time instance. We can also
+# define differences of the same matrix in different time instances and see what has changed
+a_mat_diff = Difference(a_mat, a_mat_2)
+a_dev_diff = Difference(a_dev, a_dev_2)
+a_nongpu_diff = Difference(a_nongpu, a_nongpu_2)
+
+# again we group them to a set and can show them by pl.reset(diffs_set)
+diffs_set = [a_mat_diff, a_dev_diff, a_nongpu_diff, a_mat]
diff --git a/utils/matrix_plotter/plot.py b/utils/matrix_plotter/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..16a8f1d6d0260a2bceb2bbe9a55a710e1ac32fff
--- /dev/null
+++ b/utils/matrix_plotter/plot.py
@@ -0,0 +1,446 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+class Trace:
+  """Holds all outputs collected from 1 run of the program.
+
+  This might include images of different matrices at different time instances.
+  If data has been collected only during one run, only one instance of Trace can be created
+  If, however, we are interested in comparing 2 runs (e.g. defective commit and commit which
+  we know is correct, 2 instances of Trace will be created)
+  
+  At the current moment, the only think that can be stored and viewed is the matrix stored in 
+  block cyclic distribution
+  
+  Attributes:
+    directory: path to the directory with collected data
+    name: identifier of the trace (e.g. correct_run, defective_run, etc...), just for printing purposes
+  """
+  def __init__(self, directory, name):      
+    """Initialize"""
+    self._name = name
+    self.__load_matrices(directory)
+
+  def __get_header(self, filename):
+    with open(filename, 'r') as f:
+      header = f.readline()
+      val = [int(x) for x in header.split()]
+      val.append(f.readline()[:-1])
+      print val
+    return  val
+
+  def __load_matrix(self, directory, filenames):
+    mat = []
+    for (mpi_rank, filename) in filenames:
+      #print filename, mpi_rank, prow, pcol, p_cols
+      (na, nblk, lda, localCols, my_prow, my_pcol, p_rows, p_cols, iteration, name) = self.__get_header(directory + '/' + filename)
+      if(mat == []):
+          mat = np.zeros((na, na))
+      (self._na, self._nblk, self._p_rows, self._p_cols) = (na, nblk, p_rows, p_cols)
+      prow = mpi_rank % p_rows
+      pcol = mpi_rank / p_rows
+      #print (na, nblk, lda, localCols, my_prow, my_pcol, my_nprows, my_npcols)
+      assert(my_prow == prow) 
+      assert(my_pcol == pcol) 
+        
+      loc_mat = np.loadtxt(fname = directory + '/' + filename, skiprows = 2)
+      #print "lda, localCols ", lda, localCols
+      for row_blk in range((lda-1)/nblk + 1):
+        loc_row_beg = row_blk * nblk
+        if(loc_row_beg >= lda):
+          continue
+        
+        loc_row_end = min(loc_row_beg + nblk, lda)
+              
+        glob_row_beg = (row_blk * p_rows + prow) * nblk  
+        #print "glob_row_beg = row_blk * nblk * p_rows ", glob_row_beg, row_blk, nblk, p_rows 
+        assert(glob_row_beg < na)
+        
+        glob_row_end = min(glob_row_beg + nblk, na)
+        
+        for col_blk in range((localCols-1)/nblk + 1):
+          loc_col_beg = col_blk * nblk
+          if(loc_col_beg >= localCols):
+            continue
+          
+          loc_col_end = min(loc_col_beg + nblk, localCols)
+              
+          glob_col_beg = (col_blk * p_cols + pcol) * nblk 
+          assert(glob_col_beg < na)
+          
+          glob_col_end = min(glob_col_beg + nblk, na)
+          
+          #print "local", (loc_row_beg, loc_row_end), (loc_col_beg, loc_col_end)
+          #print "global", (glob_row_beg, glob_row_end), (glob_col_beg, glob_col_end)
+          
+          mat[glob_row_beg:glob_row_end, glob_col_beg:glob_col_end] = loc_mat[loc_row_beg:loc_row_end, loc_col_beg:loc_col_end]
+    return ((name, iteration), mat)
+    
+  def __load_matrices(self, directory):
+    filenames_dict = {}
+
+    for filename in os.listdir(directory):
+      (mat_name, instant, mpi_rank) = filename.split(".")[0].split("-")
+      instant = int(instant)
+      mpi_rank = int(mpi_rank)
+
+      if(not instant in filenames_dict.keys()):
+        filenames_dict[instant] = {}
+      if(not mat_name in filenames_dict[instant].keys()):
+        filenames_dict[instant][mat_name] = []
+      filenames_dict[instant][mat_name].append((mpi_rank,filename))
+          
+    self._matrices = {}
+    self._descriptions = {}
+    for instant in filenames_dict.keys():
+      self._matrices[instant] = {}
+      for mat_name in filenames_dict[instant].keys():
+        (ids, matrix) = self.__load_matrix(directory, filenames_dict[instant][mat_name])
+        self._matrices[instant][mat_name] = matrix
+        if(not instant in self._descriptions.keys()):
+          self._descriptions[instant] = ids
+        else:
+          assert(self._descriptions[instant] == ids)
+          
+
+class Iterator:
+  """ Allows us to traverse through the history
+  
+  Contains different methods for moving forward or backword in outputed time instances
+  It allows to have different views of different matrices from (possibly) different runs
+  of the code bind to the same iterator, so that by manipulationg it (e.g. going to "next")
+  we go to next in all simultaneously active views.
+  
+  Attributes:
+    description: identifier of the iterator, to be shown in the plot
+  """
+  def __init__(self, description):
+    self._description = description
+    self._index = 0
+    self._registerd_views = []
+    
+  def __str__(self):
+    descr = self._description[self._index]
+    return "snapshot " + str(self._index) + ", iteration " + str(descr[1]) + ", " + descr[0]
+  
+  def set(self, it2):
+    """ Sets the iterator to the same state as another iterator it2 """
+    self.set_snapshot(it2._index)
+    
+  def next(self):
+    """ Moves the iterator to the next image and _updates all _attached images """
+    if(self._index < self._description.keys()[-1]):
+      self._index += 1
+      self._status_changed()
+    
+  def prev(self):
+    """ Moves the iterator to the previous image and _updates all _attached images """
+    if(self._index > 0):
+      self._index -= 1
+      self._status_changed()
+      
+  def set_snapshot(self, sn):
+    """ Sets iterator to the given snapshot sn. There can be multiple snapshots per iteration """
+    if(sn >= 0 and sn <= self._description.keys()[-1]):
+      self._index = sn
+      self._status_changed()
+      
+  def set_iteration(self, it):
+    """ Sets iterator to the given iteration """
+    for k in self._description.keys():
+      if(self._description[k][1] == it):
+        self.set_snapshot(k)
+        return
+ 
+  def next_event(self, event = None):
+    """ Sets iterator to the next event given by parameter. 
+    If not specified, moves to the next occurence of the same event it points to now"""
+    sn = self._index
+    if(event == None):
+      event = self._description[sn][0]
+    while(sn < self._description.keys()[-1]):
+      sn += 1
+      if((self._description[sn][0]).startswith(event)):
+        self.set_snapshot(sn)
+        return
+    
+  def prev_event(self, event = None):
+    """ Sets iterator to the previous event given by parameter. 
+    If not specified, moves to the previous occurence of the same event it points to now"""
+    sn = self._index
+    if(event == None):
+      event = self._description[sn][0]
+    while(sn > 0):
+      sn -= 1
+      if((self._description[sn][0]).startswith(event)):
+        self.set_snapshot(sn)
+        return   
+  
+  def get_num_iterations(self):
+    """ Returns maximal number of iteration found in the description """
+    return max([self._description[k][1] for k in self._description.keys()])
+  
+  def get_events(self):
+    """ Returns all the evetn strings. They (or their beginnings) can be used to call next(prev)_event """
+    return set([self._description[k][0] for k in self._description.keys()])
+    
+  def current(self):    
+    return self._index
+     
+  def _short_str(self):
+    descr = self._description[self._index]
+    return "sn " + str(self._index) + ", it " + str(descr[1]) + ", " + descr[0]
+  
+  def _combined_short_str(self, iter2):
+    if (self._index == iter2._index):
+      descr = self._description[self._index]
+      return "sn " + str(self._index) + ", it " + str(descr[1]) + ", " + descr[0]
+    else:
+      return "sn " + str(self._index) + "/" + str(iter2._index)
+  
+  def _register_view(self, view):
+    self._registerd_views.append(view)
+  
+  def _status_changed(self):
+    for view in self._registerd_views:
+      view._update()
+
+class View:
+  """ Abstract class representing 1 view (plot) of matrix-like data
+  It can be a matrix itself, difference of 2 matrices, etc... 
+  """
+  def __init__(self):
+    self._ax = None
+    self._show_grid = False
+  
+  def show_grid(self, show=True):
+    """ Set whether show grid between different processors """
+    self._show_grid = show
+    self._update()
+    
+  def _detach(self):
+    self._ax = None
+    self._plotter = None
+    
+  def _attach(self, plotter, ax):
+    assert(self._ax == None)
+    self._ax = ax
+    self._plotter = plotter
+    
+  def _update(self):
+    if self._ax != None:
+      self._plot(self._ax)
+      self._plotter._show()
+          
+          
+class Snapshot(View):
+  """ Represents state of 1 particular matrix from 1 particular program run (trace) at 
+  1 particular time instance (represented by state of iterator)
+  
+  It can be used to view the value of a matrix in a given time or to define difference 
+  of two snapshots
+  
+  Attributes:
+    trace: instance of Trace, determines from which run of the code we want to take data
+           (we might have saved multiple runs)
+    matrix: which matrix we want to look at (we might have saved values of more matrices)
+    iterator: determins the time instance of interest. Note that the state of Snapshot 
+              is changed by changing the state of the iterator
+  """
+  def __init__(self, trace, matrix, iterator):
+    View.__init__(self)
+    self._trace = trace
+    self._matrix = matrix
+    self._iterator = iterator    
+
+    # this is just a hack to fix the scale, it should be done better
+    # the scale is fixed to show nicely the first matrix and is kept for the rest
+    # !!! this might be very stupid for certain matrices (etc. zero at the begining )
+    self.set_limits_from_first()
+    
+    iterator._register_view(self)
+    
+  def set_limits(self, lim_lo, lim_hi):
+    """ fix limits for the plot, will be kept until changed again"""
+    self._limits = (lim_lo, lim_hi)
+    
+  def set_limits_from_first(self):
+    """fixing min and max based on first snapshot, not optimal! """
+    mat = self._trace._matrices[0][self._matrix]
+    self.set_limits(np.percentile(mat,1), np.percentile(mat,99))
+  
+  def write(self):
+    print self._matrix, ", ", it._index
+    
+  def _get_matrix(self):
+    return self._trace._matrices[self._iterator.current()][self._matrix]
+    
+  def _plot(self, ax):
+    mat = self._get_matrix().copy()
+
+    mat = np.minimum(mat, self._limits[1])
+    mat = np.maximum(mat, self._limits[0])
+    ax.clear()
+    ax.matshow(mat)
+    ax.set_title(self._trace._name + ', ' + self._matrix + ', ' + self._iterator._short_str())
+    
+    if(self._show_grid):
+      self._plot_grid(ax)
+  
+  def _plot_grid(self, ax):
+    x = 0
+    nblk = self._trace._nblk
+    na = self._trace._na
+    p_rows = self._trace._p_rows
+    p_cols = self._trace._p_cols
+    while (x < na):
+      if((x/nblk % p_rows == 0) and (x/nblk % p_cols == 0)):
+        color = 'magenta'
+      else:
+        color = 'green'
+      ax.plot([-0.5,na-0.5],[x-0.5,x-0.5], color=color, linewidth=5)
+      ax.plot([x-0.5,x-0.5],[-0.5,na-0.5], color=color, linewidth=5)
+      y = 0
+      if(x < na - 1):
+        while (y < na-1):
+          proc_row  = y/nblk % p_rows
+          proc_col  = x/nblk % p_cols
+          my_rank = proc_row + p_rows * proc_col
+          x_text = x+nblk/2-1
+          y_text = y+nblk/2+1
+          if(x_text <= na and y_text <= na):
+            ax.text(x_text, y_text , str(my_rank), fontsize = 23, color='green')
+          y+= nblk
+      x += nblk
+    
+    ax.set_xbound(-0.5,na-0.5)
+    ax.set_ybound(-0.5,na-0.5)
+    
+class Difference(View):
+  """ Compares two snapshots and shows where they differ
+  
+  Shows entries where two matrices defined by two snapshots differ
+  Intended for debugging and trying to find out where two calculations
+  start to diverge. The two snapshots can represent the same matrix in 
+  two different runs (e.g. at different git commits) or matrix on host 
+  and device, etc..
+  
+  In principle, one could define and show any operation on two matrices, 
+  the operation could be passed as callback function. At the moment, this
+  one seems to be the most reasonable
+  
+  Attributes:
+    sn1, sn2: the two snapshots to be compared
+  """
+  def __init__(self, sn1, sn2):
+    View.__init__(self)
+    self._sn1 = sn1
+    self._sn2 = sn2
+    
+    sn1._iterator._register_view(self)
+    sn2._iterator._register_view(self)
+
+  def _plot(self, ax):
+    diff = (abs(self._sn1._get_matrix() - self._sn2._get_matrix()) > 0.00001)
+    ax.clear()
+    ax.matshow(diff)
+    if(self._sn1._trace._name == self._sn2._trace._name):
+      trace_name = self._sn1._trace._name
+    else:
+      trace_name = self._sn1._trace._name + '/' +self._sn2._trace._name
+    if(self._sn1._matrix == self._sn2._matrix):
+      matrix_name = self._sn1._matrix
+    else:
+      matrix_name = self._sn1._matrix + '/' + self._sn2._matrix
+      
+    if(self._sn1._iterator._index == self._sn2._iterator._index):
+      iterator_name = self._sn1._iterator._short_str()
+    else:
+      iterator_name = self._sn1._iterator._short_str() + '/' + self._sn2._iterator._short_str()
+    
+    ax.set_title(trace_name + ', ' + matrix_name + '\n' + iterator_name)
+    if(self._show_grid):
+      self._plot_grid(ax)
+
+  def _plot_grid(self, ax):
+    self._sn1._plot_grid(ax)
+    self._sn2._plot_grid(ax)
+
+
+class Plotter:
+  """ Encapsulates the matplotlib figure and handles updating of the plot when iterators change
+  
+  It can be configured with 1 or more views by reset function. If more views are required (limited
+  to 4 at the moment), the views are automatically shown in subplots of the matplotlib plot. This 
+  can be handy when comparing more matrices at the time, especially when they share the same iterator
+  and thus one can see changes of more matrices after calling single command to the iterator.
+  
+  Attributes:
+    fig: matplotlib figure
+  """
+  def __init__(self, fig):
+    plt.ion()
+    #plt.close()    
+    self._fig = fig
+    self._views = []
+
+  def reset(self, views, title = "View"):
+    """ Defines what the Plotter should show
+    It can show 1 matrix or grid of more matrices, depending on parameter views
+    
+    Attributes:
+      views: 1 instance of View or list of up to 4 instances of View to be shown
+      title: name of the Plotter
+    """
+    if(not isinstance(views, list)):
+      self.reset([views], title)
+      return
+    
+    self._fig.clear()
+    #self._fig.tight_layout()
+    
+    self._detach_views()
+    self._views = views
+    
+    self._subplots = []
+    assert(len(views) <= 4)
+
+    if(len(views) == 1):
+      ax = self._fig.add_subplot(111)
+      ax.set_xticklabels([])
+      ax.set_yticklabels([])
+
+      self._subplots.append(ax)
+      views[0]._attach(self, ax)      
+    else:    
+      subplot_nr = 1
+      for view in views:
+        ax = self._fig.add_subplot(220 + subplot_nr)
+        ax.set_xticklabels([])
+        ax.set_yticklabels([])
+
+        self._subplots.append(ax)
+        view._attach(self, ax)      
+        subplot_nr += 1
+    
+    self._fig.suptitle(title, fontsize=20)
+    self._plot()
+    
+  def show_grid(self, show = True):
+    """ Set whether show grid between different processors, applied to all the views """
+    for view in self._views:
+      view.show_grid(show)
+
+  def _plot(self):
+    for view in self._views:
+      view._update()
+    self._fig.show()    
+    
+  def _show(self):
+    self._fig.show()
+  
+  def _detach_views(self):
+    for view in self._views:
+      view._detach()
+