void LibRaw::pseudoinverse(double (*in)[3], double (*out)[3], int size)
{
  double work[3][6], num;
  int i, j, k;

  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 6; j++)
      work[i][j] = j == i + 3;
    for (j = 0; j < 3; j++)
      for (k = 0; k < size; k++)
        work[i][j] += in[k][i] * in[k][j];
  }
  for (i = 0; i < 3; i++)
  {
    num = work[i][i];
    for (j = 0; j < 6; j++)
      if (fabs(num) > 0.00001f)
        work[i][j] /= num;
    for (k = 0; k < 3; k++)
    {
      if (k == i)
        continue;
      num = work[k][i];
      for (j = 0; j < 6; j++)
        work[k][j] -= work[i][j] * num;
    }
  }
  for (i = 0; i < size; i++)
    for (j = 0; j < 3; j++)
      for (out[i][j] = k = 0; k < 3; k++)
        out[i][j] += work[j][k + 3] * in[i][k];
}
Unknown