#include <bits/stdc++.h>
#pragma GCC optimize("Ofast")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")
using std::bitset;
using std::cin;
using std::cout;
using std::endl;
using std::make_pair;
using std::map;
using std::pair;
using std::priority_queue;
using std::queue;
using std::set;
using std::string;
using std::stringstream;
using std::vector;
class MyTimer {
 public:
  void reset() {}
  void report() {}
} myTimer;
template <typename A, typename B>
std::ostream &operator<<(std::ostream &cout, const pair<A, B> &x) {
  return cout << "(" << x.first << ", " << x.second << ")";
}
template <typename T1, typename T2>
inline bool chmin(T1 &a, const T2 &b) {
  return a > b ? a = b, true : false;
}
template <typename T1, typename T2>
inline bool chmax(T1 &a, const T2 &b) {
  return a < b ? a = b, true : false;
}
char pool[1 << 15 | 1], *it = pool + 32768;
inline int readint() {
  int a = 0;
  char c = (it >= pool + 32768
                ? (pool[fread(pool, sizeof(char), 1 << 15, stdin)] = EOF,
                   *((it = pool)++))
                : *(it++)),
       p = 0;
  while (isspace(c))
    c = (it >= pool + 32768
             ? (pool[fread(pool, sizeof(char), 1 << 15, stdin)] = EOF,
                *((it = pool)++))
             : *(it++));
  if (c == '-')
    p = 1, c = (it >= pool + 32768
                    ? (pool[fread(pool, sizeof(char), 1 << 15, stdin)] = EOF,
                       *((it = pool)++))
                    : *(it++));
  while (isdigit(c))
    a = a * 10 + c - '0',
    c = (it >= pool + 32768
             ? (pool[fread(pool, sizeof(char), 1 << 15, stdin)] = EOF,
                *((it = pool)++))
             : *(it++));
  return p ? -a : a;
}
const int maxN = 2000 + 5;
const long long P = 998244353;
int n, a[maxN][maxN];
long long f[maxN][maxN];
long long p[maxN], fac[maxN];
int C[maxN];
void clear() {
  for (int i = 0; i <= n; ++i) C[i] = 0;
}
void add(int p, int x) {
  for (; p <= n; p += p & -p) C[p] += x;
}
int sum(int p) {
  int r = 0;
  for (; p > 0; p -= p & -p) r += C[p];
  return r;
}
int main() {
  n = readint();
  for (int i = 1; i <= n; ++i)
    for (int j = 1; j <= n; ++j) a[i][j] = readint();
  f[0][0] = 1;
  f[1][0] = 0;
  f[2][0] = 1;
  for (int i = 1; i <= n; ++i) {
    if (i > 2) f[i][0] = (i - 1) * (f[i - 1][0] + f[i - 2][0]) % P;
    for (int j = 1; j <= i; ++j)
      f[i][j] = (j * f[i - 1][j - 1] + (i - j) * f[i - 1][j]) % P;
  }
  p[0] = 1;
  fac[0] = 1;
  for (int i = 1; i <= n; ++i)
    p[i] = p[i - 1] * f[n][0] % P, fac[i] = fac[i - 1] * i % P;
  long long ans = 0;
  clear();
  for (int j = 1; j <= n; ++j)
    ans = (ans + (a[1][j] - 1 - sum(a[1][j] - 1)) * fac[n - j] % P * p[n - 1]) %
          P,
    add(a[1][j], 1);
  static int posA[maxN], posB[maxN], x[maxN], y[maxN], k[maxN];
  for (int i = 2; i <= n; ++i) {
    for (int j = 1; j <= n; ++j) posA[a[i - 1][j]] = j, posB[a[i][j]] = j;
    for (int j = 1; j <= n; ++j) x[j] = y[j] = k[j] = 0;
    for (int v = 1; v <= n; ++v) {
      int lim = std::min(posA[v], posB[v]);
      for (int j = 1; j < lim; ++j) y[j] += (v < a[i][j]);
      for (int j = posA[v] + 1; j < posB[v]; ++j) x[j] += (v < a[i][j]);
      if (posA[v] <= posB[v]) k[posA[v]]++, k[posB[v] + 1]--;
    }
    for (int j = 1; j <= n; ++j) k[j] += k[j - 1];
    for (int j = 1; j <= n; ++j) {
      ans = (ans + ((x[j] == 0 ? 0 : x[j] * f[n - j][k[j] - 1]) +
                    y[j] * f[n - j][k[j]]) %
                       P * p[n - i]) %
            P;
    }
  }
  printf("%lld\n", ans % P);
  return 0;
}
