def thunk():
            in_shape = inputs[0][0].shape
            rows, cols = in_shape

            assert rows % 4 == 0

            out_shape = (rows, 4 * cols)
            
            batch_size = rows // 4
            num_features = cols

            out = outputs[0]

            # only allocate if there is no previous allocation of the right size.
            if out[0] is None or out[0].shape != out_shape:
                out[0] = cuda.CudaNdarray.zeros(out_shape)

            x_block = 16
            y_block = 16
            block = (x_block, y_block, 1)

            x_grid = int(np.ceil(float(in_shape[1]) / x_block))
            y_grid = int(np.ceil(float(in_shape[0]) / y_block))
            grid = (x_grid, y_grid, 1)

            kernel(inputs[0][0], out[0], np.intc(batch_size), np.intc(num_features), block=block, grid=grid)

    def __init__(self,
                 fc4,
                 supercell,
                 primitive,
                 mesh,
                 band_indices=None,
                 frequency_factor_to_THz=VaspToTHz,
                 is_nosym=False,
                 symprec=1e-3,
                 cutoff_frequency=1e-4,
                 log_level=False,
                 lapack_zheev_uplo='L'):
        self._fc4 = fc4
        self._supercell = supercell
        self._primitive = primitive
        self._mesh = np.intc(mesh)
        if band_indices is None:
            self._band_indices = [
                np.arange(primitive.get_number_of_atoms() * 3)]
        else:
            self._band_indices = band_indices
        self._frequency_factor_to_THz = frequency_factor_to_THz
        self._is_nosym = is_nosym
        self._symprec = symprec
        self._cutoff_frequency = cutoff_frequency
        self._log_level = log_level
        self._lapack_zheev_uplo = lapack_zheev_uplo

        self._band_indices_flatten = np.intc(
            [x for bi in self._band_indices for x in bi])

        self._frequency_shifts = None

 def x_cadd_Y_as_Y(self, alpha, x, beta, Y, result = None):
     '''       
     result[i,j] = alpha*x[j] + beta*Y[i,j]
     
     x_radd_Y_as_Y(float alpha, float* x, float beta, float* Y, float* result,
                 uint Y_col, uint Y_row)
     '''
     if result is None:
         Y_col = Y.shape[0]
         Y_row = Y.shape[1]
         self.x_cadd_Y_as_Y_kernel(np.float32(alpha), x.gpudata, \
                                   np.float32(beta), Y.gpudata, \
                                   Y.gpudata, np.intc(Y_col), np.intc(Y_row), \
                                   block = (32, 32, 1), \
                                   grid = (int(Y_row / 32) + 1, int(Y_col / 32) + 1) \
                                   )
     else:
         Y_col = Y.shape[0]
         Y_row = Y.shape[1]
         self.x_cadd_Y_as_Y_kernel(np.float32(alpha), x.gpudata, \
                                   np.float32(beta), Y.gpudata, \
                                   result.gpudata, np.intc(Y_col), np.intc(Y_row), \
                                   block = (32, 32, 1), \
                                   grid = (int(Y_row / 32) + 1, int(Y_col / 32) + 1) \
                                   )

 def _run_c(self, g_skip=None):
     import anharmonic._phono3py as phono3c
     if g_skip is None:
         g_skip = np.zeros_like(self._interaction_strength_reduced, dtype="bool")
     assert g_skip.shape == self._interaction_strength_reduced.shape
     self._set_phonon_c()
     masses = np.double(self._primitive.get_masses())
     p2s = np.intc(self._primitive.get_primitive_to_supercell_map())
     s2p = np.intc(self._primitive.get_supercell_to_primitive_map())
     atc=np.intc(self._triplet_cut_super) # int type
     atc_prim = np.intc(self._triplet_cut_prim) # int type
     phono3c.interaction(self._interaction_strength_reduced,
                         self._frequencies,
                         self._eigenvectors,
                         self._triplets_at_q_reduced.copy(),
                         self._grid_address,
                         self._mesh,
                         self._fc3,
                         atc,
                         atc_prim,
                         g_skip,
                         self._svecs,
                         self._multiplicity,
                         np.double(masses),
                         p2s,
                         s2p,
                         self._band_indices,
                         self._symmetrize_fc3_q,
                         self._cutoff_frequency,
                         self._cutoff_hfrequency,
                         self._cutoff_delta)
     phono3c.interaction_degeneracy_grid(self._interaction_strength_reduced,
                                         self._degenerates,
                                         self._triplets_at_q_reduced.astype('intc').copy(),
                                         self._band_indices.astype('intc'))

 def x_cmul_Y_as_Y(self, alpha, x, x0, Y, y0, beta, result = None):
     '''       
     result[i,j] = alpha*(x[i] + x0) * (Y[i,j] + y0) + beta
     
     x_cmul_Y_as_Y(float alpha, float* x, float x0, float* Y, float y0, float beta,
                     float* result, uint Y_col, uint Y_row)
     '''
     if result is None:
         Y_col = Y.shape[0]
         Y_row = Y.shape[1]
         self.x_cmul_Y_as_Y_kernel(np.float32(alpha), x.gpudata, np.float32(x0), \
                                   Y.gpudata, np.float32(y0), np.float32(beta), \
                                   Y.gpudata, np.intc(Y_col), np.intc(Y_row), \
                                   block = self._2d_block, \
                                   grid = self._2d_grid(Y_col, Y_row) \
                                   )
     else:
         Y_col = Y.shape[0]
         Y_row = Y.shape[1]
         self.x_cmul_Y_as_Y_kernel(np.float32(alpha), x.gpudata, np.float32(x0), \
                                   Y.gpudata, np.float32(y0), np.float32(beta), \
                                   result.gpudata, np.intc(Y_col), np.intc(Y_row), \
                                   block = self._2d_block, \
                                   grid = self._2d_grid(Y_col, Y_row) \
                                   )

 def thunk():
     
     # inputs
     A = inputs[0][0]
     B = inputs[1][0]
     
     # dimensions
     m = A.shape[0]
     n = A.shape[1]
     k = B.shape[1]
     assert n == B.shape[0] # Otherwise GEMM is impossible
     assert n%16 == 0 # Block size
     
     # output
     output_shape = (m, k)
     C = outputs[0]
     # only allocate if there is no previous allocation of the right size.
     if C[0] is None or C[0].shape != output_shape:
         C[0] = cuda.CudaNdarray.zeros(output_shape)
     
     # Launching GEMM GPU kernel            
     block_size = 16
     block = (block_size,block_size,1)
     grid = (k / block_size+1, m / block_size+1) # better too many blocks than too little
     gemm_kernel(A,B,C[0], np.intc(m), np.intc(n), np.intc(k), block= block, grid=grid)

    def set_kappa_at_s_c(self, s):
        import anharmonic._phono3py as phono3c
        kappa = np.zeros_like(self._kappa[s])
        rec_lat = np.linalg.inv(self._primitive.get_cell())
        for t, temp in enumerate(self._temperatures):
            gouterm_temp = np.zeros((self._frequencies.shape[0], self._frequencies.shape[1], 6), dtype="double")
            phono3c.phonon_gb33_multiply_dvector_gb3_dvector_gb3(gouterm_temp,
                                                                 self._b[:,t].copy(),
                                                                 self._F[s,:,t].copy(),
                                                                 np.intc(self._irr_index_mapping).copy(),
                                                                 np.intc(self._kpoint_operations[self._rot_mappings]).copy(),
                                                                 rec_lat.copy())
            kappa[:,t] = gouterm_temp * temp ** 2
            self._F[s, np.where((np.abs(self._qpoints) > self._pp._criteria).any(axis=1)), t] = 0.
            kappa[np.where((np.abs(self._qpoints) > self._pp._criteria).any(axis=1)), t] = 0.
        # l = len(np.where((np.abs(self._qpoints) > self._pp._criteria).any(axis=1))[0])
        # kappa *= np.prod(self._mesh) / np.double(np.prod(self._mesh) - l)

        kappa *= self._kappa_factor / np.prod(self._mesh)

        #modified for my own interest


        kappa_max = kappa.sum(axis=(0,2)).max(axis=-1)
        rkappa = np.sum(np.abs(kappa - self._kappa[s]), axis=(0, 2)) # over qpoints and nbands
        for i in range(6):
            self._rkappa[s, :, i] = rkappa[:, i] /  kappa_max
        self._kappa[s] = kappa

 def setNeighborPair(self, s1, s2, w):
     """Create an edge (s1, s2) with weight w.
     w should be of type int or float (python primitive type).
     s1 should be smaller than s2."""
     if not (0 <= s1 < s2 < self.numSites):
         raise IndexOutOfBoundError()
     _cgco.gcoSetNeighborPair(self.handle, np.intc(s1), np.intc(s2), self._convertPairwiseTerm(w))

    def route_flow(self, receiver, dem='topographic__elevation'):
        #main
        self._dem = self._grid['node'][dem]
        """
        if receiver==None:
            self._flow_receiver = self._flow_dirs_d8(self._dem)
        else:
            self._flow_receiver = receiver
        """
        self._flow_receiver = receiver
        #(self._flow_receiver, ss) = grid_flow_directions(self._grid, self._dem)

        method = 'python'
        if method=='cython':
            from flow_direction_over_flat_cython import adjust_flow_direction
            self._flow_receiver = adjust_flow_direction(self._dem, np.intc(self._flow_receiver),
                                                        np.intc(self._boundary), np.intc(self._close_boundary),
                                                        np.intc(self._open_boundary), np.intc(self._neighbors))
        else:
            flat_mask, labels = self._resolve_flats()
            self._flow_receiver = self._flow_dirs_over_flat_d8(flat_mask, labels)


        #a, q, s = flow_accum_bw.flow_accumulation(self._flow_receiver, self._open_boundary, node_cell_area=self._grid.forced_cell_areas)

        #self._grid['node']['flow_receiver'] = self._flow_receiver

        return self._flow_receiver

def get_triplets_reciprocal_mesh_at_q(fixed_grid_number,
                                      mesh,
                                      rotations,
                                      is_time_reversal=True,
                                      is_return_map=False,
                                      is_return_rot_map=False):

    weights = np.zeros(np.prod(mesh), dtype='intc')
    third_q = np.zeros(np.prod(mesh), dtype='intc')
    mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
    mapping = np.zeros(np.prod(mesh), dtype='intc')
    rot_mapping = np.zeros(np.prod(mesh), dtype='intc')

    spg.triplets_reciprocal_mesh_at_q(weights,
                                      mesh_points,
                                      third_q,
                                      mapping,
                                      rot_mapping,
                                      fixed_grid_number,
                                      np.intc(mesh).copy(),
                                      is_time_reversal * 1,
                                      np.intc(rotations).copy())
    assert len(mapping[np.unique(mapping)]) == len(weights[np.nonzero(weights)]), \
        "At grid %d, number of irreducible mapping: %d is not equal to the number of irreducible triplets%d"%\
            (fixed_grid_number, len(mapping[np.unique(mapping)]), len(weights[np.nonzero(weights)]))
    if not is_return_map and not is_return_rot_map:
        return weights, third_q, mesh_points
    elif not is_return_rot_map:
        return weights, third_q,mesh_points, mapping
    else:
        return weights, third_q,mesh_points, mapping, rot_mapping

    def __init__(self,
                 fc4,
                 supercell,
                 primitive,
                 mesh,
                 temperatures=None,
                 band_indices=None,
                 frequency_factor_to_THz=VaspToTHz,
                 is_nosym=False,
                 symprec=1e-3,
                 cutoff_frequency=1e-4,
                 log_level=False,
                 lapack_zheev_uplo='L'):
        self._fc4 = fc4
        self._supercell = supercell
        self._primitive = primitive
        self._masses = np.double(self._primitive.get_masses())
        self._mesh = np.intc(mesh)
        if temperatures is None:
            self._temperatures = np.double([0])
        else:
            self._temperatures = np.double(temperatures)
        num_band = primitive.get_number_of_atoms() * 3
        if band_indices is None:
            self._band_indices = np.arange(num_band, dtype='intc')
        else:
            self._band_indices = np.intc(band_indices)
        self._frequency_factor_to_THz = frequency_factor_to_THz
        self._is_nosym = is_nosym
        self._symprec = symprec
        self._cutoff_frequency = cutoff_frequency
        self._log_level = log_level
        self._lapack_zheev_uplo = lapack_zheev_uplo

        symmetry = Symmetry(primitive, symprec=symprec)
        self._point_group_operations = symmetry.get_pointgroup_operations()

        self._grid_address = None
        self._bz_map = None
        self._set_grid_address()
        
        self._grid_point = None
        self._quartets_at_q = None
        self._weights_at_q = None

        self._phonon_done = None
        self._frequencies = None
        self._eigenvectors = None
        self._dm = None
        self._nac_q_direction = None

        self._frequency_shifts = None
        
        # Unit to THz of Delta
        self._unit_conversion = (EV / Angstrom ** 4 / AMU ** 2
                                 / (2 * np.pi * THz) ** 2
                                 * Hbar * EV / (2 * np.pi * THz) / 8
                                 / np.prod(self._mesh))

        self._allocate_phonon()

 def X_router_Y_add_O(self, alpha, X, Y, beta, O, result = None):
     '''
     xi, yi is vector
     X, Y is matrix
     
     X=[x1, x2, x3]
     Y=[y1, y2, y3]
     X_router_Y_add_O = a * (x1 outer y1 + x2 outer y2 + x3 outer y3) + b * O
     
     X_router_Y_add_O(float alpha, float* X, float* Y, float beta, float* O, float* result, uint r_col,
     uint r_row, uint XY_col)      
     '''
     if result is None:
         O_col = O.shape[0]
         O_row = O.shape[1]
         XY_col = X.shape[0]
         self.X_router_Y_add_O_kernel(np.float32(alpha), X.gpudata, Y.gpudata, \
                                      np.float32(beta), O.gpudata, \
                                      O.gpudata, np.intc(O_col), np.intc(O_row), np.intc(XY_col), \
                                      block = (32, 32, 1), \
                                      grid = (int(O_row / 32) + 1, int(O_col / 32) + 1) \
                                      )
     else:
         O_col = O.shape[0]
         O_row = O.shape[1]
         XY_col = X.shape[0]
         self.X_router_Y_add_O_kernel(np.float32(alpha), X.gpudata, Y.gpudata, \
                                      np.float32(beta), O.gpudata, \
                                      result.gpudata, np.intc(O_col), np.intc(O_row), np.intc(XY_col), \
                                      block = (32, 32, 1), \
                                      grid = (int(O_row / 32) + 1, int(O_col / 32) + 1) \
                                      )

    def _calculate_fc4_normal_c(self):
        import anharmonic._phono4py as phono4c
        svecs, multiplicity = get_smallest_vectors(self._supercell,
                                                   self._primitive,
                                                   self._symprec)
        p2s = np.intc(self._primitive.get_primitive_to_supercell_map())
        s2p = np.intc(self._primitive.get_supercell_to_primitive_map())
        gp = self._grid_point
        self._set_phonon_c([gp])
        self._set_phonon_c(self._quartets_at_q)

        phono4c.fc4_normal_for_frequency_shift(
            self._fc4_normal,
            self._frequencies,
            self._eigenvectors,
            gp,
            self._quartets_at_q,
            self._grid_address,
            self._mesh,
            np.double(self._fc4),
            svecs,
            multiplicity,
            self._masses,
            p2s,
            s2p,
            self._band_indices,
            self._cutoff_frequency)

def get_csr_matrix(A):
    '''get csr matrix from dolfin without copying data

    cf. http://code.google.com/p/pyamg/source/browse/branches/2.0.x/Examples/DolfinFormulation/demo.py
    '''
    (row,col,data) = A.data()
    return csr_matrix( (data,intc(col),intc(row)), shape=(A.size(0),A.size(1)) )

def get_reduced_triplets_permute_sym(triplets,
                                     mesh,
                                     first_mapping,
                                     first_rotation,
                                     second_mapping):
    mesh = np.array(mesh)
    triplet_numbers = np.array([len(tri) for tri in triplets], dtype="intc")
    grid_points = np.array([triplet[0][0] for triplet in triplets], dtype="intc")
    triplets_all = np.vstack(triplets)
    triplets_mapping = np.arange(len(triplets_all)).astype("intc")
    sequences = np.array([[0,1,2]] * len(triplets_all), dtype="byte")

    num_irred_triplets =  spg.reduce_triplets_permute_sym(triplets_mapping,
                                           sequences,
                                           triplets_all.astype("intc"),
                                           np.intc(grid_points).copy(),
                                           np.intc(triplet_numbers).copy(),
                                           mesh.astype("intc"),
                                           np.intc(first_mapping).copy(),
                                           np.intc(first_rotation).copy(),
                                           np.intc(second_mapping).copy())
    assert len(np.unique(triplets_mapping)) == num_irred_triplets
    unique_triplets, indices = np.unique(triplets_mapping, return_inverse=True)
    triplets_map = []
    triplet_sequence = []
    num_triplets = 0
    for i, triplets_at_q in enumerate(triplets):
        triplets_map.append(indices[num_triplets:num_triplets+len(triplets_at_q)])
        triplet_sequence.append(sequences[num_triplets:num_triplets+len(triplets_at_q)])
        num_triplets += len(triplets_at_q)
    return unique_triplets,triplets_map, triplet_sequence

    def calculate_gh_at_sigma_and_temp(self):
        import anharmonic._phono3py as phono3c
        if self._is_precondition:
            out = self._collision_out[self._isigma, :, self._itemp]
            out_reverse = np.where(self._frequencies>self._cutoff_frequency, 1 / out, 0)
            self._z[self._isigma, :, self._itemp] = self._r[self._isigma, :, self._itemp] * out_reverse[..., np.newaxis]
        else:
            self._z[self._isigma, :, self._itemp] = self._r[self._isigma, :, self._itemp]

        self._z[:, np.where(np.any(np.abs(self._qpoints) > self._pp._criteria, axis=1))] = 0
        zr0 = np.zeros(3, dtype="double")
        phono3c.phonon_3_multiply_dvector_gb3_dvector_gb3(zr0,
                                                        self._z_prev[self._isigma,:,self._itemp].copy(),
                                                        self._r_prev[self._isigma,:,self._itemp].copy(),
                                                        np.intc(self._irr_index_mapping).copy(),
                                                        np.intc(self._kpoint_operations[self._rot_mappings]),
                                                        np.double(np.linalg.inv(self._primitive.get_cell())).copy())
        #Flexibly preconditioned CG: r(i+1)-r(i) instead of r(i+1)
        r = self._r[self._isigma,:,self._itemp] - self._r_prev[self._isigma,:,self._itemp]
        # r = self._r[self._isigma,:,self._itemp]
        zr1 = np.zeros(3, dtype="double")
        phono3c.phonon_3_multiply_dvector_gb3_dvector_gb3(zr1,
                                                        self._z[self._isigma,:,self._itemp].copy(),
                                                        r.copy(),
                                                        np.intc(self._irr_index_mapping).copy(),
                                                        np.intc(self._kpoint_operations[self._rot_mappings]),
                                                        np.double(np.linalg.inv(self._primitive.get_cell())).copy())
        zr1_over_zr0 = np.where(np.abs(zr0>0), zr1/zr0, 0)
        self._p[self._isigma,:,self._itemp] = self._z[self._isigma,:,self._itemp] +\
                                              zr1_over_zr0 * self._p_prev[self._isigma,:,self._itemp]

        self._p[:, np.where(np.any(np.abs(self._qpoints) > self._pp._criteria, axis=1))] = 0

def get_reduced_pairs_permute_sym(pairs,
                                 mesh,
                                 first_mapping,
                                 first_rotation,
                                 second_mapping):
    mesh = np.array(mesh)
    pair_numbers = np.array([len(pair) for pair in pairs], dtype="intc")
    grid_points = np.array([pair[0][0] for pair in pairs], dtype="intc")
    pairs_all = np.vstack(pairs)
    pairs_mapping = np.arange(len(pairs_all)).astype("intc")
    sequences = np.array([[0,1]] * len(pairs_all), dtype="byte")

    num_irred_pairs =  spg.reduce_pairs_permute_sym(pairs_mapping,
                                           sequences,
                                           pairs_all.astype("intc"),
                                           np.intc(grid_points).copy(),
                                           np.intc(pair_numbers).copy(),
                                           mesh.astype("intc"),
                                           np.intc(first_mapping).copy(),
                                           np.intc(first_rotation).copy(),
                                           np.intc(second_mapping).copy())
    assert len(np.unique(pairs_mapping)) == num_irred_pairs
    unique_pairs, indices = np.unique(pairs_mapping, return_inverse=True)
    pairs_map = []
    pair_sequence = []
    num_pairs = 0
    for i, pairs_at_q in enumerate(pairs):
        pairs_map.append(indices[num_pairs:num_pairs+len(pairs_at_q)])
        pair_sequence.append(sequences[num_pairs:num_pairs+len(pairs_at_q)])
        num_pairs += len(pairs_at_q)
    return unique_pairs,pairs_map, pair_sequence

def get_mappings(mesh,
                 rotations,
                 is_shift=np.zeros(3, dtype='intc'),
                 is_time_reversal=True,
                 qpoints=np.double([])):
    """
    Return k-point map to the irreducible k-points and k-point grid points .

    The symmetry is searched from the input rotation matrices in real space.

    is_shift=[0, 0, 0] gives Gamma center mesh and the values 1 give
    half mesh distance shifts.
    """

    mapping = np.zeros(np.prod(mesh), dtype='intc')
    rot_mapping = np.zeros(np.prod(mesh), dtype="intc")
    mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
    qpoints = np.double(qpoints).copy()
    if qpoints.shape == (3,):
        qpoints = np.double([qpoints])
    spg.stabilized_reciprocal_mesh(mesh_points,
                                   mapping,
                                   rot_mapping,
                                   np.intc(mesh).copy(),
                                   np.intc(is_shift),
                                   is_time_reversal * 1,
                                   np.intc(rotations).copy(),
                                   np.double(qpoints))

    return mapping,  rot_mapping

def get_kpoint_group(mesh,  point_group_operations, qpoints=[[0.,0.,0.]], is_time_reversal=True):
    reverse_kpt_group = np.zeros((len(point_group_operations) * 2, 3, 3), dtype="intc")
    num_kpt_rots = spg.kpointgroup(reverse_kpt_group,
                                    np.intc(point_group_operations).copy(),
                                    np.intc(mesh).copy(),
                                    np.double(qpoints).copy(),
                                    is_time_reversal)
    return reverse_kpt_group[:num_kpt_rots]

def to_scipy_csr(mat, dtype=None, imagify=False):
    (row,col,data) = mat.data()   # get sparse data
    col = intc(col)
    row = intc(row)
    n = mat.size(0)
    if imagify: data = data*1j
    Asp = csr_matrix( (data,col,row), shape=(n,n), dtype=dtype)
    return Asp