13

The singular value decomposition for a symmetric matrix $A=A^{T}$ is one and the same as its canonical eigendecomposition (i.e. with an orthonormal matrix-of-eigenvectors), while the same thing for a nonsymmetric matrix $M=U \Sigma V^T$ is just the canonical eigenvalue decomposition for the symmetric matrix $$ H=\begin{bmatrix}0 & M\\ M^{T} & 0 \end{...


9

There's no reason to append a row of 1's. You should just perform a rank-revealing QR factorization (like with routine SGEQP3) on $A^T$, and the last column of $Q$ should be in the nullspace. This has the added advantage that the relative magnitude of the last element on the diagonal of $R$ gives you some idea of how singular the solution is. Even better ...


9

This is usually caused by trying to use a threaded MKL combined with MPI, resulting in over-subscription. Either explicitly configure PETSc to use non-threaded MKL or add MKL_NUM_THREADS=1 to your environment.


7

Prior answers to this question have covered most of the salient points, but I want to add one comment with respect to this: does MKL have the upper hand for some tasks? The MKL team is in a unique position to know about future Intel instruction sets and their implementations in specific processors. Furthermore, they have access to proprietary processor ...


5

The various Fortran standards allow a lot of compiler dependent behaviour in terms of function binary interfaces when being called with "complicated" data types such as Fortran90 style arrays and complex numbers. This means calling code compiled with one compiler from another one is not guaranteed to do what you expect, and can lead to grabbing the wrong bit ...


4

I believe that MKL has the threaded parallel and serial functions in one unified library. You can try setting OMP_NUM_THREADS or MKL_NUM_THREADS to a range of values and see how the performance varies. Setting either to 1 will give you the serial behavior.


3

There is some good information here: https://software.intel.com/en-us/mkl-linux-developer-guide-calling-intel-mkl-functions-from-multi-threaded-applications set MKL_NUM_THREADS=n call mkl_set_num_threads(n)


3

IPS and IPC are generally specified "per core". That's because processor makers often vary how many cores of a particular kind they pack on the same processor, so it doesn't really make sense to specify these per-processor, whereas the core is always the same in those cases -- the type of core is generally described by the "generation" of the processor. I ...


3

From the old document: Intel® Math Kernel Library FFT to DFTI Wrappers, A314775-001US: All transforms require additional memory to store the transform coefficients. When performing multiple FFTs of the same dimension, the table of coefficients should be created only once and then used on all the FFTs afterwards. Using the same table rather than creating ...


3

Nothing stops you from decomposing the problem up yourself and feeding the relevant partitioned data into MKL sequentially, or even in parallel. It will work as long as you avoid data races, but you may experience performance penalties unless you are very careful about how you do it. The reason it's discouraged to combine OpenMP code with MKL is that OpenMP ...


2

Note that Armadillo's cx_double (which is used by cx_vec) is just an alias for std::complex<double>, so in fact you're converting between std::complex<double> and MKL_Complex16. The only requirement for the reinterpret_cast to behave as expected is the compatibility of memory layouts of the involved structures -- in this case they are compatible. ...


2

Looking at the MKL documentation, I see that MKL has only CG and GMRES solvers, along with ILU preconditioners. It may be that the GMRES solver will be adequate for your purposes.


2

A matrix of size 15M x 15M is likely too big for a (sparse) direct solver on a single machine -- it is going to take too much time and memory. If you wanted to use a direct solver, you could try parallel sparse direct solvers such as MUMPS or SuperLU-dist, both of which are conveniently called via PETSc. There is also the option of using iterative solvers. ...


2

Being a big believer in parallel computation (the chips are not getting any faster), the short answer is yes. You should look into MPI (Message Passing Interface), the defecto standard for passing data between different machines when it comes to cluster computing. See this link https://en.wikipedia.org/wiki/Message_Passing_Interface for an intro to MPI. ...


1

https://software.intel.com/en-us/mkl-developer-reference-c-mkl-sparse-syrk This routine is specifically designed for your problem. It will output the upper half of the resulting matrix, which is often preferable.


1

Although the question has a great answer, here's a rule of thumb for small singular values, with a plot. If a singular value is nonzero but very small, you should define its reciprocal to be zero, since its apparent value is probably an artifact of roundoff error, not a meaningful number. A plausible answer to the question "how small is small ?" ...


1

Here is the Eigen example you requested. It follows the approach in Christian Clason's comment and also does the same computation using an explicit inverse for comparison. #include <iostream> using std::cout; using std::endl; #include <ctime> #include <Eigen/Core> #include <Eigen/LU> void compareInvAndFactor() { typedef Eigen::...


1

There is source code of this function available here: http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html As you can see by following the function calls in the code, it computes the Cholesky factorization of $A$ with dpotrf, followed by dpotrs, which solves $AX=B$ by solving two triangular systems (dtrsm). The implementations of ...


1

Computing the condition number of a matrix is an expensive operation (for example, involving computing the largest and smallest magnitude eigenvalues, or solving several(!) linear systems), so doing this by default in a call to a linear solver would be very wasteful. It makes much more sense to just try to solve the system, and in case the solver breaks down ...


1

If this is a large band matrix, I assume it is sparse and you could well use the sparse solvers in MKL (such as Pardiso or iterative such as conjugate gradients) http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-78889273-7E77-426A-9B5E-23A7C2378D78.htm


1

Use the MKL link line advisor to determine what library that corresponds to. https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor If you choose "Single Dynamic Library", it uses "mkl_rt.lib" and it looks like threading is active.


1

Summarizing the answer to the question. zgetrf (LU factorization of a complex matrix $A$) returns info=i>0 if $U_{ii}=0$; thus in your case $U_{44}=0$ during the attempted LU-decomposition by MKL. That happens (excluding very exotic cases) when your original matrix $A$ is singular or numerically singular. If zgetrf returned an error there is no point ...


Only top voted, non community-wiki answers of a minimum length are eligible