LedHead
  • Member for 4 years, 11 months
  • Last seen more than 1 year ago
  • Columbus, OH
How to start using LAPACK in c++?
Accepted answer
24 votes

I'm going to disagree with some of the other answers and say that I believe that figuring out how to use LAPACK is important in the field of scientific computing. However, there is a large learning ...

View answer
Writing a non-square linear system in standard form $A\cdot{x}=b$
Accepted answer
8 votes

The matrix math in this paper is terribly hard to follow, but here's what equation (19) should look like $\left( \begin{array}{c|c|c|c} c_{11} \mathbf{I} - \mathbf{A}^T & c_{12} \mathbf{I} &...

View answer
Sparse matrix inversion
Accepted answer
7 votes

For a matrix that small, you're probably not going to do better than using dense methods. I wrote up a quick test in C++ for an 18x18 matrix with your sparse structure and randomly generated values ...

View answer
Parametrized spline - oscilating second derivative
Accepted answer
7 votes

Your choice of parameterization is creating problems. Instead of spanning one in $t$ between points, span an amount proportional to the line segment between the two points in $(x,y)$ space. I've ...

View answer
Floating Point error when computing Binomial Distribution Probability
Accepted answer
5 votes

This is why gammaln exists logB = gammaln(m+n+1) - gammaln(m+1) - gammaln(n+1) + m*log(p) + n*log(q);

View answer
Right-preconditioning and fixed point linear iterations
Accepted answer
5 votes

Starting with $AM^{-1} \mathbf{y} = \mathbf{b}$, where $\mathbf{y} \equiv M \mathbf{x}$, we can manipulate to $\mathbf{y} - (I-AM^{-1})\mathbf{y} = \mathbf{b}$ Replacing one instance of $\mathbf{...

View answer
Transform from linear index of a packed triangular matrix to dense indices
Accepted answer
4 votes

Given the linear index $f$, you can invert $f(i,j) = i+\frac{j(j+1)}{2}$ by first calculating the $j$-index via $j = \left\lfloor \frac{-1+\sqrt{1+8f}}{2} \right\rfloor$ where $\lfloor \cdot \...

View answer
Working with large mesh files
4 votes

My answer is primarily opinion-based, given my experience. In my work, I haven't (yet) dealt with meshes quite as large as what you're describing. However, I've seen large enough meshes to hint that ...

View answer
Plot a multi-variable equation in (discontinuous) regions
Accepted answer
3 votes

First, you originally wrote that your equation is $\cos(K)=5 \text{ sinc}(5.12\sqrt{E})- \cos(5.12\sqrt{E})$, but you clearly meant $\cos(K)=5 \text{ sinc}(5.12\sqrt{E}) + \cos(5.12\sqrt{E})$ Second, ...

View answer
What exactly is Hardware Acceleration?
3 votes

Hardware acceleration means any code run on specialized hardware, as opposed to software run on general purpose CPUs such as standard x86 processors on your PC. I suppose the term is inherently ...

View answer
A good, simple book/resource on Parallel Programming in C++ for scientific computing
3 votes

Here's a couple of resources for MPI for C language. mpitutorial.com A User's Guide to MPI by Peter Pacheco. This is accessible from the page that @BillGreene referenced in his comment, but this is ...

View answer
Parallelizing evolutionary algorithm on PBS
2 votes

I think what you're trying to do is totally reasonable. I don't think there's any problem running your ea_loop script on the head node provided it's doing mostly control level operations like ...

View answer
What are some algorithms to calculate the width of an arbitrary polygon when a bounding box approximation is inaccurate
2 votes

I'm sure there are better solutions than this, but since no one else has answered to this point, I'll throw out a this-is-what-I'd-do answer. Triangulate the polygon If your polygon doesn't have too ...

View answer
Rewriting matrix multiplication
2 votes

Solution slightly modified from here W2 = N * bsxfun(@times, N, G).'; This works for N of size m x n, and G of size 1 x n.

View answer
numpy.outer without flatten
2 votes

I think numpy.tensordot does what you need. import numpy as np N=2 M=3 L=4 x=np.arange(N*M).reshape(N,M) y=np.arange(L) z=np.tensordot(x,y,axes=0) print('x=',x) print('y=',y) print('z=',z) x= [[0 ...

View answer
Naming scientific output files
2 votes

This question might get closed for being opinion-based, but I think it's an important question ... so here's my opinion. As far as I know, there's no real industry best practice for this kind of ...

View answer
How to broadcast my matrix exponentiation
2 votes

Since $A$ is only a 2x2 matrix, you could just hard-code the matrix exponentiation. Note that if we can write $A=UDU^{-1}$, where $D$ is a diagonal matrix of eigenvalues of $A$ and the columns of $U$ ...

View answer
problem in interface operator
1 votes

Here's the flow of logic, hopefully in a slightly more readable form: We write the transmission conditions with as-yet-unknown linear operators $\mathcal{S}_1$ and $\mathcal{S}_2$ which operate on ...

View answer
Is there an efficient way to form this block matrix with numpy or scipy?
1 votes

If you need to explicitly construct the entire matrix, then Stefano M's answer is your best bet. If, however, you don't really need the whole matrix, but just need to be able to perform a matrix-...

View answer
Visually appealing ways to plot singular vector fields with matplotlib or other foss tools
1 votes

I think it makes sense to present your data in logarithmic scale for both color and vector length. Here's a code snippet based on your example. from matplotlib.colors import LogNorm from pylab import ...

View answer
Electromagnetism FEM (FEniCS) interpolation - leakage effect
0 votes

Your problem is using project() in FEniCS. Here's the FEniCS documentation that discusses why you might want to use the project operator. Note that in that example, the exact flux is continuous, while ...

View answer
How to use LAPACK function (DGELSY) in Fortran
Accepted answer
0 votes

You're most of the way there. On input, a is the matrix A and b is the right-hand side vector or matrix b. When the subroutine finishes, a has been completely changed (it now holds the QR ...

View answer
Specifying mesh spacing for DFT in numpy
0 votes

I'd love to say that I understand exactly all of the scale factors and shifts that I did below, but I mostly played around with factors until things matched :) I would wait until I've thought it all ...

View answer
Optimization techniques for expensive multi-variable functions
0 votes

You're asking for a method, but let me point you to a specific software package - SNOBFIT. I don't think it can handle discrete variables, but it does handle noisy cost function evaluations, which is ...

View answer
Exploiting patterns in matrix for efficient matrix-vector multiplication
0 votes

Disclaimer: I have no idea if this will actually speed up your calculation, as it adds quite a bit of computational overhead. Since it looks like your matrix isn't very sparse, it's hard to imagine ...

View answer
find grid points inside the parallelogram defined by an origin and two vectors
0 votes

I think you're just looking for a straight-forward way of mapping from your original space, call it $(x,y)$ space to your $(a,b)$ space, where $(x,y)=(i\Delta x, j\Delta y)$. Once you're in $(a,b)$ ...

View answer
Underdetermined/overdetermined polynomial interpolation
Accepted answer
0 votes

I'm assuming you're doing your interpolation something like the following way: Given a set of $(x,y,z)$ coordinates $(x_i,y_i,z_i)$ and corresponding value $f_i$ for $i=1,2,...,N$, find the ...

View answer