/**
  * created     : 30.10.2024 12:36:13
  * author      : wzj33300
  */

#include <bits/stdc++.h>
using namespace std;

#ifdef DEBUG
#include <algo/debug.h>
#else
#define debug(...) 42
#define assert(...) 42
#endif

using ll = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using db = long double;  // or double, if TL is tight
using str = string;      // yay python!

// pairs
using pi = pair<int, int>;
using pl = pair<ll, ll>;
using pd = pair<db, db>;
#define mp make_pair
#define fi first
#define se second

// ^ lol this makes everything look weird but I'll try it
template <class T>
using vc = vector<T>;
template <class T, size_t SZ>
using AR = array<T, SZ>;
using vi = vc<int>;
using vb = vc<bool>;
using vl = vc<ll>;
using vd = vc<db>;
using vs = vc<str>;
using vpi = vc<pi>;
using vpl = vc<pl>;
using vpd = vc<pd>;

// vectors
// oops size(x), rbegin(x), rend(x) need C++17
#define sz(x) int((x).size())
#define bg(x) begin(x)
#define all(x) bg(x), end(x)
#define rall(x) x.rbegin(), x.rend()
#define sor(x) sort(all(x))
#define rsz resize
#define ins insert
#define pb push_back
#define eb emplace_back
#define ft front()
#define bk back()

#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep1(i, n) for (int i = 1; i < (n); ++i)
#define rep1n(i, n) for (int i = 1; i <= (n); ++i)
#define repr(i, n) for (int i = (n) - 1; i >= 0; --i)

#define rep_subset(t, s) \
  for (ll t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s)))

#define lb lower_bound
#define ub upper_bound
template <class T>
int lwb(vc<T> &a, const T &b) {
  return int(lb(all(a), b) - bg(a));
}
template <class T>
int upb(vc<T> &a, const T &b) {
  return int(ub(all(a), b) - bg(a));
}
// const int MOD = 998244353;  // 1e9+7;
const int Big = 1e9;  // not too close to INT_MAX
const ll BIG = 1e18;  // not too close to LLONG_MAX
const db PI = acos((db)-1);
const int dx[4]{1, 0, -1, 0}, dy[4]{0, 1, 0, -1};  // for every grid problem!!
mt19937 rng((uint32_t)chrono::steady_clock::now().time_since_epoch().count());
template <class T>
using pqg = priority_queue<T, vector<T>, greater<T>>;

int pct(int x) { return __builtin_popcount(x); }
int pct(u32 x) { return __builtin_popcount(x); }
int pct(ll x) { return __builtin_popcountll(x); }
int pct(u64 x) { return __builtin_popcountll(x); }
int popcnt_mod_2(int x) { return __builtin_parity(x); }
int popcnt_mod_2(u32 x) { return __builtin_parity(x); }
int popcnt_mod_2(ll x) { return __builtin_parityll(x); }
int popcnt_mod_2(u64 x) { return __builtin_parityll(x); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2)
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(u32 x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 0, 2)
int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(u32 x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
int lowbit(u64 x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }

template <typename T>
T floor(T a, T b) {
  return a / b - (a % b && (a ^ b) < 0);
}
template <typename T>
T ceil(T x, T y) {
  return floor(x + y - 1, y);
}
template <typename T>
T bmod(T x, T y) {
  return x - y * floor(x, y);
}
template <typename T>
pair<T, T> divmod(T x, T y) {
  T q = floor(x, y);
  return {q, x - q * y};
}

template <class T>
bool ckmin(T &a, const T &b) {
  return b < a ? a = b, 1 : 0;
}  // set a = min(a,b)
template <class T>
bool ckmax(T &a, const T &b) {
  return a < b ? a = b, 1 : 0;
}  // set a = max(a,b)

template <class T, class U>
T fstTrue(T lo, T hi, U f) {
  ++hi;
  assert(lo <= hi);  // assuming f is increasing
  while (lo < hi) {  // find first index such that f is true
    T mid = lo + (hi - lo) / 2;
    f(mid) ? hi = mid : lo = mid + 1;
  }
  return lo;
}
template <class T, class U>
T lstTrue(T lo, T hi, U f) {
  --lo;
  assert(lo <= hi);  // assuming f is decreasing
  while (lo < hi) {  // find first index such that f is true
    T mid = lo + (hi - lo + 1) / 2;
    f(mid) ? lo = mid : hi = mid - 1;
  }
  return lo;
}
int mod1, mod2;
int pw1[400 * 400 + 1], pw2[400 * 400 + 1];
struct Hash {
  int x1 = 0, x2 = 0;
  friend Hash &operator+=(Hash &me, int x) {
    me.x1 += pw1[x];
    me.x1 %= mod1;
    me.x2 += pw2[x];
    me.x2 %= mod2;
    return me;
  }
  friend Hash operator+(const Hash &x, const Hash &y) {
    return Hash{(x.x1 + y.x1) % mod1, (x.x2 + y.x2) % mod2};
  }
  friend bool operator!=(const Hash &x, const Hash &y) {
    return x.x1 != y.x1 || x.x2 != y.x2;
  }
};

const int N = 400 * 400 * 20 * 2,
          mod = 1e9 + 7;

int ls[N], rs[N], cnt, val[N], pw[400 * 400 + 1], c[N];
Hash sum[N];

void update(int &rt, int l, int r, int x) {
  ++cnt;
  ls[cnt] = ls[rt], rs[cnt] = rs[rt], sum[cnt] = sum[rt], val[cnt] = val[rt], c[cnt] = c[rt];
  rt = cnt;
  c[rt]++;
  val[rt] += pw[x];
  val[rt] %= mod;
  sum[rt] += x;
  if (l == r) {
    return;
  }
  int mid = l + r >> 1;
  if (x <= mid)
    update(ls[rt], l, mid, x);
  else
    update(rs[rt], mid + 1, r, x);
}

bool cmp(int x, int y, int l, int r) {
  if (l == r) return c[x] > c[y];
  int mid = l + r >> 1;
  if (sum[rs[x]] != sum[rs[y]])
    return cmp(rs[x], rs[y], mid + 1, r);
  else
    return cmp(ls[x], ls[y], l, mid);
}

using P = pair<int, int>;

P LS(P x) {
  return P{ls[x.fi], ls[x.se]};
}
P RS(P x) {
  return P{rs[x.fi], rs[x.se]};
}
Hash SUM(P x) {
  return sum[x.fi] + sum[x.se];
}
int VAL(P x) {
  return (val[x.fi] + val[x.se]) % mod;
}
int C(P x) {
  return c[x.fi] + c[x.se];
}

bool cmp(P x, P y, int l, int r) {
  if (l == r) return C(x) > C(y);
  int mid = l + r >> 1;
  if (SUM(RS(x)) != SUM(RS(y)))
    return cmp(RS(x), RS(y), mid + 1, r);
  else
    return cmp(LS(x), LS(y), l, mid);
}

void _sol() {
  int n, m;
  cin >> n >> m;
  pw[0] = 1;
  for (int i = 1; i <= n * n; i++) {
    pw[i] = 1ll * pw[i - 1] * n * n % mod;
  }
  vc<vi> a(n, vi(n));
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
      cin >> a[i][j];
  vc<vc<int>> f(n, vc<int>(n));
  cnt = 0;
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++) {
      if (i == 0 && j == 0) {
      } else if (i == 0)
        f[i][j] = f[i][j - 1];
      else if (j == 0)
        f[i][j] = f[i - 1][j];
      else
        f[i][j] = cmp(f[i - 1][j], f[i][j - 1], 1, n * n) ? f[i - 1][j] : f[i][j - 1];
      update(f[i][j], 1, n * n, a[i][j]);
      // debug(val[f[i][j]]);
    }
  vc<vc<int>> g(n, vc<int>(n));
  vc<vc<P>> pre(n, vc<P>(n)), suf(pre);
  repr(i, n) repr(j, n) {
    if (i == n - 1 && j == n - 1) {
    } else if (i == n - 1)
      g[i][j] = g[i][j + 1];
    else if (j == n - 1)
      g[i][j] = g[i + 1][j];
    else
      g[i][j] = cmp(g[i + 1][j], g[i][j + 1], 1, n * n) ? g[i + 1][j] : g[i][j + 1];
    pre[i][j] = suf[i][j] = P{f[i][j], g[i][j]};
    update(g[i][j], 1, n * n, a[i][j]);
  }
  // debug(f, g);
  for (int i = 0; i < n; i++) {
    for (int j = 1; j < n; j++) {
      pre[i][j] = cmp(pre[i][j], pre[i][j - 1], 1, n * n) ? pre[i][j] : pre[i][j - 1];
    }
    for (int j = n - 2; j >= 0; j--) {
      suf[i][j] = cmp(suf[i][j], suf[i][j + 1], 1, n * n) ? suf[i][j] : suf[i][j + 1];
    }
  }
  while (m--) {
    int x1, x2, y1, y2;
    cin >> x1 >> x2 >> y1 >> y2;
    --x1, --x2, --y1, --y2;
    P ans{0, 0};
    // (x2, y1)
    if (x2 < n - 1 && y1 > 0) {
      ans = cmp(ans, pre[x2 + 1][y1 - 1], 1, n * n) ? ans : pre[x2 + 1][y1 - 1];
    }
    // (x1, y2)
    if (x1 > 0 && y2 < n - 1) {
      ans = cmp(ans, suf[x1 - 1][y2 + 1], 1, n * n) ? ans : suf[x1 - 1][y2 + 1];
    }
    cout << VAL(ans) << '\n';
  }
}

// signed main() {
int main() {
  // freopen("ex_d1.in", "r", stdin);
  // freopen(".out","w",stdout);
  ios::sync_with_stdio(false);
  cin.tie(0);
  mod1 = 998244353;
  mod2 = 1e9 + 7;
  for (int i = 0; i <= 400 * 400; i++) {
    pw1[i] = rng() % mod1;
    pw2[i] = rng() % mod2;
  }
  int t;
  cin >> t;
  while (t--) {
    _sol();
  }
  return 0;
}