Open WalterLu3 opened 7 months ago
Thanks for the detailed reproducible code. I was able to run it here. I'll upload my files as a PR, if you like.
The matrix can be factorized, but when used as is stands, the matrix takes about 151 seconds on my 20-core desktop. It suffers a lot of fillin.
I tried pulling the matrix into MATLAB (which also uses umfpack), and I noticed something interesting. The matrix has 530,379 entries but only 448,277 nonzeros. That is, 82,102 of the entries are exactly zero. Normally that's not a problem, but I do treat those as legitimate entries (it can happen that subsequent matrices need those entries).
But when I tried it in MATLAB, the factorization took just 2 seconds, and L+U has just 19.6 million entries. MATLAB drops zeros by design.
If I add the zeros with values equal to epsilon (1e-16), the factorization takes 258 seconds in MATLAB (an older version of UMFPACK perhaps explains that difference, or a different BLAS library), and the factorization L+U has 1 billion entries. I think the epsilon values I added are causing trouble there. I'm not factorizing the same matrix in MATLAB as in C, in this case.
I don't see as many entries in L+U from the C code (where zeros are kept but equal to zero), but I still see a lot: about 90 million, and the time is 151 seconds.
So I would consider dropping the explicit zeros from A, if you can do that.
I'll try that and let you know what happens. Thanks.
Here's my output with your original matrix:
that's my output, with some added diagnostics, with this program:
loadUMFPACK_long_version.c.txt
and a build script: CMakeLists.txt
my MATLAB script: try_matlab.m.txt
and the output:
and its output: diary.txt
My desktop has 256GB of RAM, so it was able to handle this matrix.
Hi, after dropping the zeros. It is able to factorize with just 32-bit interface on my 16GB RAM machine!
#include <stdio.h> //perror
#include "umfpack.h"
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
/*********************************************************************************
* failed_matrix.dat is a binary that stores our instances.
*
* failed_matrix.dat contains the matrix information in the
* following order :
*
* 1. int size (the size of the square matrix. i.e the number of rows and columns)
* 2. int nnz (number of non-zero element in our data matrix)
* 3. int *Ai (An array of length 'nnz' to store the row indices
* 4. int *Ap (An array of length 'size + 1')
* 5. double *Ax (An array of length 'nnz' to store matrix data)
*
*
* So it stores everything based on UMFPACK's form.
* *******************************************************************************/
/*********************************************************************************
* This file does the following steps
* 1. Read in the matrix that fails to be factorized.
* 2. Dynamically load the routines for factorization. (you need to input your
* UMFPACK_PATH to link to the dynamic libraray)
* 3. Try to factorize it using int32 interface.
* 4. Failed with Out-Of-Memory Error
********************************************************************************/
int main(int argc, char **argv){
FILE *fptr;
int status;
int size;
int nnz;
int *Ai;
int *Ap;
double *Ax;
int *Ai2;
int *Ap2;
double *Ax2;
double *null = (double *) NULL ;
void *Symbolic, *Numeric ;
/* read in matrix */
/* open file*/
fptr = fopen("failed_matrix.dat","r");
if (fptr == NULL){
perror("file open failed");
}
/* read in all matrix information */
fread(&size, sizeof(int), 1, fptr);
fread(&nnz, sizeof(int), 1, fptr);
Ai = (int *)malloc(nnz*sizeof(int));
Ap = (int *)malloc((size + 1)*sizeof(int));
Ax = (double *)malloc(nnz*sizeof(double));
fread(Ai, sizeof(int), nnz, fptr);
fread(Ap, sizeof(int), size + 1, fptr);
fread(Ax, sizeof(double), nnz, fptr);
/* finish reading matrix*/
Ai2 = (int *)malloc(nnz*sizeof(int));
Ap2 = (int *)malloc((size + 1)*sizeof(int));
Ax2 = (double *)malloc(nnz*sizeof(double));
int number_nnz = 0;
int col_nnz;
int colStart, colEnd;
int start;
/* drop zeros in the matrix */
Ap2[0] = 0;
for (int i = 0; i <=size; i++){
col_nnz = 0;
start = Ap[i];
colStart = Ap[i];
colEnd = Ap[i+1];
while(colStart < colEnd){
if (Ax[colStart] != 0){
Ai2[number_nnz] = Ai[colStart];
Ax2[number_nnz] = Ax[colStart];
number_nnz++;
col_nnz++;
}
colStart ++;
}
Ap2[i + 1] = Ap2[i] + col_nnz;
}
printf("%d,%d\n",Ap2[size + 1],number_nnz);
printf("zeros : %d\n", number_nnz);
/* start factorization step*/
printf("run symbolic step\n");
status = umfpack_di_symbolic(size, size, Ap2, Ai2, Ax2, &Symbolic, null, null);
if (status != 0){
perror("symbolic");
}
printf("run numeric step\n");
status = umfpack_di_numeric(Ap2, Ai2, Ax2, Symbolic, &Numeric, null, null);
if (status != 0){
printf("status code : %d\n", status);
}
free(Ai);
free(Ap);
free(Ax);
return 0;
}
Maybe HiGHS has done something for the zero elements in the matrix, I will probably look more into that when I have time. Thanks for the detailed explanation!
Glad to help!
Can you tell me more about this matrix? It could make for an interesting problem to add to my collection, as a counter-example.
Sure. I was doing work on adding new basis packages to the PATH solver, which is a specialized solver for the mixed complementarity problem. The basis packages are used to calculate the basic feasible direction for pivoting since the PATH solver adopts a method similar to the simplex method. When I was testing all the available packages in PATH with some problem sets, a problem failed (a dynamic equilibrium problem). The failed matrix I extracted corresponds to the basic feasible direction, $B^{-1}A_j$ of the pivot where things collapse.
The reason why it has zeros in the data is exactly what you said - it can happen that subsequent matrices need those entries (we indeed reuse the matrix structure). Essentially, we are doing things similar to extracting the columns of Jacobian of a nonlinear map to form a basis. When we are identifying the number of non-zero elements, we are identifying the entries that are "structurally nonzero." So an entry that has value $x_1x_2$ will be identified as non-zero, but the value may be zero because it depends on the given value of $x_1$ and $x_2$.
Describe the bug I encountered an instance of matrix where UMFPACK 32-bit routine is unable to factorize it because of out-of-memory error. I then tried to use the UMFPACK 64-bit routine but it get stuck at the _numeric stage. I am able to factorize this matrix using other basis routines such as the ones provided by HiGHS. Not sure if it’s a bug, but I hope reporting this is helpful for you.
To Reproduce The following are the steps I do to reproduce the error. The binary file that contains the matrix is failed_matrix.dat. All the files mentioned here can be found in: https://github.com/WalterLu3/umfpack_edge_case
Step 1. Reproducing Out-Of-Memory error with UMFPACK 32-bit interface.
I compiled this by dynamically linking to libumfpack.dylib.
It gives me an out-of-memory error. According to the UMFPACK user guide, I proceeded to use UMPACK 64 routines.
Step 2. Reproducing error with UMFPACK 64-bit never ending. I also compiled the following by dynamically linking to libumfpack.dylib.
The numeric steps seem to never stops.
Expected behavior I expected the numeric step would end for the UMFPACK 64bit interface.
Desktop (please complete the following information):
Thanks.