code: def diag(self, reordered = True):
          
          sncolptr = self.symb.sncolptr
          snptr = self.symb.snptr
          snode = self.symb.snode
          blkptr = self.symb.blkptr
          
          D = matrix(0.0,(self.symb.n,1))
          for k in range(self.symb.Nsn):
              nn = snptr[k+1]-snptr[k]
              w = sncolptr[k+1]-sncolptr[k]
              for i in range(nn): D[snode[snptr[k]+i]] = self.blkval[blkptr[k]+i*w+i]
          if reordered: return D
          else: return D[self.symb.ip]
docs: Returns a vector with the diagonal elements of the matrix.