#define _USE_MATH_DEFINES
#define eps 1.0e-12
#define mz 1.0e-15
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <algorithm>
#include <assert.h>
#include <time.h>
#include <ctime>
using namespace std;
void qralgShift(vector<double>* matrix, int dim, vector<double>* resids);
//qr factorization
extern “C” void dgeqrf_(const int *M, const int *N, double *A, const int *lda,
double *tau, double *work, int *lwork, int *info);

//return new matrix in matrix
//return residuals in resids
//to get number of qr factorizations, look at how many resids there are
//shift is Wilkinson shift
void qralgShift(vector* matrix, int dim, vector* resids)
{
int i,j,k,l;

//stuff for lapack qr function
int M = dim; //number of rows
int N = dim; //number of cols
double *A;
int lda = M; //the leading dimension of the array A
double *tau; //min(M,N)
double *work;
int lwork = N;
int info;

A = (double*)malloc(M*N*sizeof(double));
tau = (double*)malloc(M*sizeof(double));
work = (double*)malloc(lwork*sizeof(double));

vectorR(dim*dim);
vectorQ(dim*dim);
vectornewQ(dim*dim);
vectorHQ(dim*dim);
vectorv(dim);
vector vvt(dim*dim);
vectorH(dim*dim);
double mu, sigma, SIG;
double resid;

double a, b, c,d,e;

do
{
//determine shift
// lower corner is:
// a(m-1) b(m-1)
// b(m-1) a(m)

//(a(m-1)-a(m))/2
sigma = 0.5*((*matrix)[(dim-2)*dim + (dim-2)] – (*matrix)[(dim-1)*dim + (dim-1)]);
SIG = -1.0;
if(sigma > mz)
SIG = 1.0;

// a(m) – (( sign(sigma) * b(m-1)^2) / (|sigma| + sqrt( sigma^2 + b(m-1)^2) )
mu = (*matrix)[(dim-1)*dim + (dim-1)] – ((SIG * (*matrix)[(dim-1)*dim + (dim-2)] * (*matrix)[(dim-1)*dim + (dim-2)]) / (fabs(sigma) + sqrt((sigma*sigma) + ((*matrix)[(dim-1)*dim + (dim-2)] * (*matrix)[(dim-1)*dim + (dim-2)]))));

//subtract shift from diagonals (A-muI)
for(j=0; j<dim; j++)//rows
{
for(k=0; k<dim; k++)//cols
{
if(j==k)
{
(*matrix)[j*dim+k] -= mu;
}
}
}

//printMatrix(dim, dim, &(*matrix));

//copy matrix into A, col-major; this is so I can give the matrix to the qr alg
for(i=0; i<dim; i++) //row matrix, col A
{
for(j=0; j= n);
//**************************************************************
R.assign(dim*dim,0.0); //initialize R to zeros
//copy A diagonal and above in R (this is R)
for(i=0; i<dim; i++) //row matrix, col A
{
for(j=0; j=i)
R[i*dim+j] = A[j*dim+i];
}
}
//I have R!

//**************************************************************
//Q: the elements below the diagonal,
//with the array TAU, represent the orthogonal matrix Q as a
//product of min(m,n) elementary reflectors
//The matrix Q is represented as a product of elementary reflectors
//Q = H(1) H(2) . . . H(k), where k = min(m,n).
//Each H(i) has the form
//H(i) = I – tau * v * vT
//where tau is a real scalar, and v is a real vector with
//v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
//and tau in TAU(i).
//**************************************************************
Q.assign(dim*dim, 0.0);
//copy A below diagonal in Q
for(i=0; i<dim; i++) //row matrix, col A
{
for(j=0; j<dim; j++) //col matrix, row A
{
if(j==i)
Q[i*dim+j] = 1.0;
if(j<i)
Q[i*dim+j] = A[j*dim+i];
}
}

//now cols of Q are my vectors v
for(i=0; i<dim; i++)
{
//get ith v
for(j=0; j<dim; j++)
v[j] = Q[j*dim+i];

//find vvT
for(j=0; j<dim; j++)
{
for(k=0; k<dim; k++)
{
vvt[j*dim+k] = v[j]*v[k];
}
}

//find I – tau[i]* vvT
H.assign(dim*dim, 0.0);
for(j=0; j<dim; j++)//rows
{
for(k=0; k<dim; k++)//cols
{
if(j==k)
{
H[j*dim+k] = 1.0 – tau[i]*vvt[j*dim+k];
}
else
{
H[j*dim+k] = -tau[i]*vvt[j*dim+k];
}

}
}

if(i==0)
{
//copy H into newQ
for(j=0; j<dim*dim; j++)
newQ[j] = H[j];
}
else
{
//multiply newQ and H
for(j=0; j<dim; j++)
{
for(k=0; k<dim; k++)
{
HQ[j*dim+k] = 0.0;
for(l=0; l<dim; l++)
{
HQ[j*dim+k] += newQ[j*dim+l] * H[l*dim+k];
}
}
}
//copy HQ into newQ
for(j=0; j<dim*dim; j++)
newQ[j] = HQ[j];
}

}
//now I have Q!

//newQ*R should give A – it does
//R*newQ is the next qr matrix
for(j=0; j<dim; j++)
{
for(k=0; k<dim; k++)
{
(*matrix)[j*dim+k] = 0.0;
for(l=0; l<dim; l++)
{
(*matrix)[j*dim+k] += R[j*dim+l] * newQ[l*dim+k];
}
}
}

//add shift to diagonals (A+muI)
for(j=0; j<dim; j++)//rows
{
for(k=0; k eps);

free (A);
}

A nice concise explanation of BVPs and the shooting method + the algorithm for the same.

Shooting Method PDF

Scenario: your code is computing partial derivatives (gradients) of some function F with respect to one or more variables. You are getting numbers out, but you are not sure if they are correct.
Solution: use finite-difference approximations to check partial derivatives

1. Take the variable with respect to which you are computing partial derivative (if you have several, then pick one) – call it q.
2. Pick a constant (call it ε). It should be a value from e-4 to e-8.
3. Change the value of q to q+ε (perturb up). Compute value of your function F using this new q_up. Call it F_up.
4. Change the value of q to q-ε (perturb down). Compute value of your function F using q_down. Call it F_down.
5. Now the central finite difference approximation of the partial derivative for the variable q is (F_up-F_down) / (2.0*epsilon).
   Compare this with the number your code is producing  (partial derivative of q w.r.t. F). The approximation should be within 10% of your code’s value, unless the derivatives are close to zero (machine zero).
   If your numbers are very different,  try changing the value of ε (if you used e-5, try using e-4 or e-6).
   If you cannot get your numbers to agree, then make yourself some coffee and start looking for bugs.

References: Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization