LU() [1/4]

M2_arrayintOrNull Lapack::LU ( const DMatCC * A, DMatCC * L, DMatCC * U )

static

Definition at line 1117 of file lapack.cpp.

{
  int rows = static_cast<int>(A->numRows());
  int cols = static_cast<int>(A->numColumns());
  int info;
  int min = (rows <= cols) ? rows : cols;
  M2_arrayint result = M2_makearrayint(rows);
 
  L->resize(rows, min);
  U->resize(min, cols);
 
  if (min == 0)
    {
      if (rows > 0)
        for (int i = 0; i < rows; i++) result->array[i] = i;
      return result;
    }
 
  int* perm = new int[min];
  std::vector<double> copyA = make_lapack_array(*A);
 
  zgetrf_(&rows, &cols, copyA.data(), &rows, perm, &info);
 
  if (info < 0)
    {
      ERROR("argument passed to zgetrf had an illegal value");
      result = nullptr;
    }
  else
    {
      fill_lower_and_upper(copyA, *L, *U);
 
      for (int i = 0; i < rows; i++) result->array[i] = i;
      for (int i = 0; i < min; i++)
        {
          int thisloc = perm[i] - 1;
          int tmp = result->array[thisloc];
          result->array[thisloc] = result->array[i];
          result->array[i] = tmp;
        }
 
// TODO: What is this block?  Remove it?
#if 0
      /* set the permutation array */
      for (int row=1; row<=min; row++) {
        int targ = row;
        for (int i=1; i<=min; i++) {
          if (i == targ)
            targ = perm[i-1];
          else if (perm[i-1] == targ)
            targ = i;
        }
        result->array[row-1] = targ-1;
      }
      for (int i=min; i<rows; i++)
        result->array[i] = i;
#endif
    }
 
  delete[] perm;
  return result;
}

References ERROR, fill_lower_and_upper(), M2_makearrayint(), make_lapack_array(), DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), DMat< ACoeffRing >::resize(), result(), and zgetrf_().