#include<stdio.h>

/* solves least squares Ax=b with
 * A     1 2
 *       3 4
 *       5 6
 *  b   3
 *      7
 *      11
 *
 *      obviously, the solution is 
 *      1
 *      1
 *
 *
 *      */
void  dqrdc_(double *A,int *ldx,int *n,int *p,double *qraux,int *jpvt,double *work, int *job);
void  dqrsl_(double *A,int *ldA,int *n,int *p,double *qraux,double *b,double *,double *,double *ans,double *,double *,int *job,int *info);


int main(int argc, char ** argv)
{
	//first dimension the arrays - note fortran uses column major
	//rather than row major
	//note also that FORTRAN calls by reference so always give addresses/pointers!!
#define ACOL 2
#define AROW 3

double A[ACOL][AROW]= { {1, 3, 5}, {2,4,6}};
double b[AROW]={3, 7, 11};

//now set up auxiliary vars for fortran routines


int jpvt[ACOL];//for pivots - but we won't pivot!!
int ldA=AROW,n=AROW,p=ACOL,job=0; //no pivot
double* work=NULL; //need to allocate if were to pivot
double qraux[ACOL];
int info;
double ans[ACOL];
double qty[AROW];
  //first to factorise
 dqrdc_(A,&ldA,&n,&p,qraux,jpvt,work,&job);

 job = 100; //for dqrsl - flag to compute least squares
  //now to solve
 dqrsl_(A,&ldA,&n,&p,qraux,b,NULL,qty,ans,NULL,NULL,&job,&info);
 
	 //print out answer
 printf("solution is (%g,%g)\n",ans[0],ans[1]);
}