<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Danyang :</div><div class="gmail_quote"><br></div><div class="gmail_quote">Adding '-mat_superlu_dist_fact SamePattern' fixed the problem. Below is how I figured it out.</div><div class="gmail_quote"><br></div><div class="gmail_quote"><div class="gmail_quote">1. Reading ex52f.F, I see '-superlu_default' = '-pc_factor_mat_solver_package superlu_dist', the later enables runtime options for other packages. I use superlu_dist-4.2 and superlu-4.1 for the tests below.</div><div class="gmail_quote"><br></div><div class="gmail_quote">2. Use the Matrix 168 to setup KSP solver and factorization, all packages, petsc, superlu_dist and mumps give same correct results:</div><div class="gmail_quote"><br></div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_168.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package petsc</div><div class="gmail_quote"> -->loac matrix a</div><div class="gmail_quote"> -->load rhs b</div><div class="gmail_quote"> size l,m,n,mm       90000       90000       90000       90000</div><div class="gmail_quote">Norm of error  7.7308E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">..</div><div class="gmail_quote"> -->Test for matrix          172</div><div class="gmail_quote">Norm of error  3.8461E-11 iterations     1</div><div class="gmail_quote"><br></div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_168.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package superlu_dist</div><div class="gmail_quote">Norm of error  9.4073E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">...</div><div class="gmail_quote"> -->Test for matrix          172</div><div class="gmail_quote">Norm of error  3.8187E-11 iterations     1</div><div class="gmail_quote"><br></div><div class="gmail_quote"><div class="gmail_quote">3. Use superlu, I get</div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_168.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package superlu</div><div class="gmail_quote">Norm of error  1.0191E-06 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">...</div><div class="gmail_quote"> -->Test for matrix          172</div><div class="gmail_quote">Norm of error  9.7858E-07 iterations     1</div><div class="gmail_quote"><br></div><div class="gmail_quote">Replacing default DiagPivotThresh: 1. to 0.0, I get same solutions as other packages:</div><div class="gmail_quote"><br></div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_168.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package superlu -mat_superlu_diagpivotthresh 0.0</div><div class="gmail_quote"><br></div><div class="gmail_quote">Norm of error  8.3614E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">...</div><div class="gmail_quote"> -->Test for matrix          172</div><div class="gmail_quote">Norm of error  3.7098E-11 iterations     1</div><div class="gmail_quote"><br></div><div class="gmail_quote">4.</div><div class="gmail_quote">using '-mat_view ascii::ascii_info', I found that a_flow_check_1.bin and a_flow_check_168.bin seem have same structure:</div><div class="gmail_quote"><br></div><div class="gmail_quote"> -->loac matrix a</div><div class="gmail_quote">Mat Object: 1 MPI processes</div><div class="gmail_quote">  type: seqaij</div><div class="gmail_quote">  rows=90000, cols=90000</div><div class="gmail_quote">  total: nonzeros=895600, allocated nonzeros=895600</div><div class="gmail_quote">  total number of mallocs used during MatSetValues calls =0</div><div class="gmail_quote">    using I-node routines: found 45000 nodes, limit used is 5</div><div class="gmail_quote"><br></div><div class="gmail_quote"><div class="gmail_quote">5.</div><div class="gmail_quote">Using a_flow_check_1.bin, I am able to reproduce the error you reported: all packages give correct results except superlu_dist:</div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_1.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package superlu_dist</div><div class="gmail_quote">Norm of error  2.5970E-12 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">Norm of error  1.3936E-01 iterations    34</div><div class="gmail_quote"> -->Test for matrix          169</div><div class="gmail_quote"><br></div><div class="gmail_quote">I guess the error might come from reuse of matrix factor. Replacing default</div><div class="gmail_quote">-mat_superlu_dist_fact <SamePattern_SameRowPerm> with</div><div class="gmail_quote">-mat_superlu_dist_fact SamePattern, I get</div><div class="gmail_quote"><br></div><div class="gmail_quote">./ex52f -f0 matrix_and_rhs_bin/a_flow_check_1.bin -rhs matrix_and_rhs_bin/b_flow_check_168.bin -loop_matrices flow_check -loop_folder matrix_and_rhs_bin -pc_type lu -pc_factor_mat_solver_package superlu_dist -mat_superlu_dist_fact SamePattern</div><div class="gmail_quote"><br></div><div class="gmail_quote">Norm of error  2.5970E-12 iterations     1</div><div class="gmail_quote"> -->Test for matrix          168</div><div class="gmail_quote">Norm of error  9.4073E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          169</div><div class="gmail_quote">Norm of error  6.4303E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          170</div><div class="gmail_quote">Norm of error  7.4327E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          171</div><div class="gmail_quote">Norm of error  5.4162E-11 iterations     1</div><div class="gmail_quote"> -->Test for matrix          172</div><div class="gmail_quote">Norm of error  3.4440E-11 iterations     1</div><div class="gmail_quote"> --> End of test, bye</div><div class="gmail_quote"><br></div><div class="gmail_quote">Sherry may tell you why SamePattern_SameRowPerm cause the difference here.</div><div class="gmail_quote">Best on the above experiments, I would set following as default</div><div class="gmail_quote">'-mat_superlu_diagpivotthresh 0.0' in petsc/superlu interface.</div><div class="gmail_quote">'-mat_superlu_dist_fact SamePattern' in petsc/superlu_dist interface.</div><div class="gmail_quote"><br></div><div class="gmail_quote">Hong</div></div></div></div><div class="gmail_quote"><br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
  
    
  
  <div bgcolor="#FFFFFF" text="#000000">
    Hi Hong,<br>
    <br>
    I did more test today and finally found that the solution accuracy
    depends on the initial (first) matrix quality. I modified the
    ex52f.F to do the test. There are 6 matrices and right-hand-side
    vectors. All these matrices and rhs are from my reactive transport
    simulation. Results will be quite different depending on which one
    you use to do factorization. Results will also be different if you
    run with different options. My code is similar to the First or the
    Second test below. When the matrix is well conditioned, it works
    fine. But if the initial matrix is well conditioned, it likely to
    crash when the matrix become ill-conditioned. Since most of my case
    are well conditioned so I didn't detect the problem before. This
    case is a special one. <br>
    <br>
    <br>
    How can I avoid this problem? Shall I redo factorization? Can PETSc
    automatically detect this prolbem or is there any option available
    to do this?<br>
    <br>
    All the data and test code (modified ex52f) can be found via the
    dropbox link below. <br>
    <u><br>
    </u><u><a href="https://www.dropbox.com/s/4al1a60creogd8m/petsc-superlu-test.tar.gz?dl=0" target="_blank">https://www.dropbox.com/s/4al1a60creogd8m/petsc-superlu-test.tar.gz?dl=0</a></u><br>
    <br>
    <br>
    Summary of my test is shown below. <br>
    <br>
    First, use the Matrix 1 to setup KSP solver and factorization, then
    solve 168 to 172<br>
    <br>
    mpiexec.hydra -n 1 ./ex52f -f0
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/a_flow_check_1.bin
    -rhs
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/b_flow_check_1.bin
    -loop_matrices flow_check -loop_folder
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin -pc_type lu
    -pc_factor_mat_solver_package superlu_dist<br>
    <br>
    Norm of error  3.8815E-11 iterations     1<br>
     -->Test for matrix          168<br>
    Norm of error  4.2307E-01 iterations    32<br>
     -->Test for matrix          169<br>
    Norm of error  3.0528E-01 iterations    32<br>
     -->Test for matrix          170<br>
    Norm of error  3.1177E-01 iterations    32<br>
     -->Test for matrix          171<br>
    Norm of error  3.2793E-01 iterations    32<br>
     -->Test for matrix          172<br>
    Norm of error  3.1251E-01 iterations    31<br>
    <br>
    Second, use the Matrix 1 to setup KSP solver and factorization using
    the implemented SuperLU relative codes. I thought this will generate
    the same results as the First test, but it actually not.<br>
    <br>
    mpiexec.hydra -n 1 ./ex52f -f0
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/a_flow_check_1.bin
    -rhs
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/b_flow_check_1.bin
    -loop_matrices flow_check -loop_folder
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin
    -superlu_default<br>
    <br>
    Norm of error  2.2632E-12 iterations     1<br>
     -->Test for matrix          168<br>
    Norm of error  1.0817E+04 iterations     1<br>
     -->Test for matrix          169<br>
    Norm of error  1.0786E+04 iterations     1<br>
     -->Test for matrix          170<br>
    Norm of error  1.0792E+04 iterations     1<br>
     -->Test for matrix          171<br>
    Norm of error  1.0792E+04 iterations     1<br>
     -->Test for matrix          172<br>
    Norm of error  1.0792E+04 iterations     1<br>
    <br>
    <br>
    Third, use the Matrix 168 to setup KSP solver and factorization,
    then solve 168 to 172<br>
    <br>
    mpiexec.hydra -n 1 ./ex52f -f0
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/a_flow_check_168.bin
    -rhs
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/b_flow_check_168.bin
    -loop_matrices flow_check -loop_folder
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin -pc_type lu
    -pc_factor_mat_solver_package superlu_dist<br>
    <br>
    Norm of error  9.5528E-10 iterations     1<br>
     -->Test for matrix          168<br>
    Norm of error  9.4945E-10 iterations     1<br>
     -->Test for matrix          169<br>
    Norm of error  6.4279E-10 iterations     1<br>
     -->Test for matrix          170<br>
    Norm of error  7.4633E-10 iterations     1<br>
     -->Test for matrix          171<br>
    Norm of error  7.4863E-10 iterations     1<br>
     -->Test for matrix          172<br>
    Norm of error  8.9701E-10 iterations     1<br>
    <br>
    Fourth, use the Matrix 168 to setup KSP solver and factorization
    using the implemented SuperLU relative codes. I thought this will
    generate the same results as the Third test, but it actually not.<br>
    <br>
    mpiexec.hydra -n 1 ./ex52f -f0
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/a_flow_check_168.bin
    -rhs
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin/b_flow_check_168.bin
    -loop_matrices flow_check -loop_folder
    /home/dsu/work/petsc-superlu-test/matrix_and_rhs_bin
    -superlu_default<br>
    <br>
    Norm of error  3.7017E-11 iterations     1<br>
     -->Test for matrix          168<br>
    Norm of error  3.6420E-11 iterations     1<br>
     -->Test for matrix          169<br>
    Norm of error  3.7184E-11 iterations     1<br>
     -->Test for matrix          170<br>
    Norm of error  3.6847E-11 iterations     1<br>
     -->Test for matrix          171<br>
    Norm of error  3.7883E-11 iterations     1<br>
     -->Test for matrix          172<br>
    Norm of error  3.8805E-11 iterations     1<br>
    <br>
    Thanks very much,<br>
    <br>
    Danyang<br>
    <br>
    <div>On 15-12-03 01:59 PM, Hong wrote:<br>
    </div>
    <blockquote type="cite">
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote">Danyang :</div>
          <div class="gmail_quote">Further testing a_flow_check_168.bin,</div>
          <div class="gmail_quote">
            <div class="gmail_quote">./ex10 -f0
              /Users/Hong/Downloads/matrix_and_rhs_bin/a_flow_check_168.bin
              -rhs
              /Users/Hong/Downloads/matrix_and_rhs_bin/x_flow_check_168.bin
              -pc_type lu -pc_factor_mat_solver_package superlu
              -ksp_monitor_true_residual -mat_superlu_conditionnumber</div>
            <div class="gmail_quote">  Recip. condition number =
              1.610480e-12</div>
            <div class="gmail_quote">  0 KSP preconditioned resid norm
              6.873340313547e+09 true resid norm 7.295020990196e+03
              ||r(i)||/||b|| 1.000000000000e+00</div>
            <div class="gmail_quote">  1 KSP preconditioned resid norm
              2.051833296449e-02 true resid norm 2.976859070118e-02
              ||r(i)||/||b|| 4.080672384793e-06</div>
            <div class="gmail_quote">Number of iterations =   1</div>
            <div class="gmail_quote">Residual norm 0.0297686</div>
            <div class="gmail_quote"><br>
            </div>
            <div class="gmail_quote">condition number of this matrix =
              1/1.610480e-12 = 1.e+12,</div>
            <div class="gmail_quote">i.e., this matrix is
              ill-conditioned. </div>
            <div class="gmail_quote"><br>
            </div>
            <div class="gmail_quote">Hong</div>
            <div class="gmail_quote"><br>
            </div>
            <div class="gmail_quote"><br>
            </div>
            <blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
              <div text="#000000" bgcolor="#FFFFFF"> Hi Hong,<br>
                <br>
                The binary format of matrix, rhs and solution can be
                downloaded via the link below.<br>
                <br>
                <a href="https://www.dropbox.com/s/cl3gfi0s0kjlktf/matrix_and_rhs_bin.tar.gz?dl=0" target="_blank">https://www.dropbox.com/s/cl3gfi0s0kjlktf/matrix_and_rhs_bin.tar.gz?dl=0</a><br>
                <br>
                Thanks,<br>
                <br>
                Danyang
                <div>
                  <div><br>
                    <br>
                    On 15-12-03 10:50 AM, Hong wrote:<br>
                    <blockquote type="cite">
                      <div dir="ltr">
                        <div class="gmail_extra">
                          <div class="gmail_quote">Danyang:<br>
                            <blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
                              <div text="#000000" bgcolor="#FFFFFF"><br>
                                <br>
                                To my surprising, solutions from SuperLU
                                at timestep 29 seems not correct for the
                                first 4 Newton iterations, but the
                                solutions from iteration solver and
                                MUMPS are correct. <br>
                                <br>
                                Please find all the matrices, rhs and
                                solutions at timestep 29 via the link
                                below. The data is a bit large so that I
                                just share it through Dropbox. A piece
                                of matlab code to read these data and
                                then computer the norm has also been
                                attached. <br>
                                <u><a href="https://www.dropbox.com/s/rr8ueysgflmxs7h/results-check.tar.gz?dl=0" target="_blank">https://www.dropbox.com/s/rr8ueysgflmxs7h/results-check.tar.gz?dl=0</a></u></div>
                            </blockquote>
                            <div><br>
                            </div>
                            <div>Can you send us matrix in petsc binary
                              format?</div>
                            <div><br>
                            </div>
                            <div>e.g., call MatView(M,
                              PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD))</div>
                            <div>or '-ksp_view_mat binary'</div>
                            <div><br>
                            </div>
                            <div>Hong</div>
                            <blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
                              <div text="#000000" bgcolor="#FFFFFF"><br>
                                <br>
                                Below is a summary of the norm from the
                                three solvers at timestep 29, newton
                                iteration 1 to 5.<br>
                                <br>
                                Timestep 29<br>
                                Norm of residual seq 1.661321e-09,
                                superlu 1.657103e+04, mumps 3.731225e-11
                                <br>
                                Norm of residual seq 1.753079e-09,
                                superlu 6.675467e+02, mumps 1.509919e-13
                                <br>
                                Norm of residual seq 4.914971e-10,
                                superlu 1.236362e-01, mumps 2.139303e-17
                                <br>
                                Norm of residual seq 3.532769e-10,
                                superlu 1.304670e-04, mumps 5.387000e-20
                                <br>
                                Norm of residual seq 3.885629e-10,
                                superlu 2.754876e-07, mumps 4.108675e-21
                                <br>
                                <br>
                                Would anybody please check if SuperLU
                                can solve these matrices? Another
                                possibility is that something is wrong
                                in my own code. But so far, I cannot
                                find any problem in my code since the
                                same code works fine if I using
                                iterative solver or direct solver MUMPS.
                                But for other cases I have tested,  all
                                these solvers work fine.<br>
                                <br>
                                Please let me know if I did not write
                                down the problem clearly.<br>
                                <br>
                                Thanks,<br>
                                <br>
                                Danyang<br>
                                <br>
                                <br>
                                <br>
                              </div>
                            </blockquote>
                          </div>
                          <br>
                        </div>
                      </div>
                    </blockquote>
                    <br>
                  </div>
                </div>
              </div>
            </blockquote>
          </div>
          <br>
        </div>
      </div>
    </blockquote>
    <br>
  </div>

</blockquote></div><br></div></div>