Sparse, Banded Matrix Solution using Mac OS X cLAPACK

gistfile1.m

- (BOOL) simulateWorkspace
/*"
**	This method will take the existing workspace with all the objects placed
**	on it and simulate it for the potential at each simulation grid point.
**	This needs to be done before you can get any values out of the workspace,
**	but that's pretty obvious if you think about it.
"*/
{
	BOOL			error = NO;

	// first, make sure we're set up for this
	if (!error) {
		if (([self getRowCount] <= 1) || ([self getColCount] <= 1) ||
			([self getEpsilonR] == nil) || ([self getRho] == nil) ||
			([self getVoltage] == nil)) {
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - this workspace is not yet set up properly for a simulation. You need to set reasonable simulation node counts as well as initializing this class for the simulation. Please make sure you cann one of the -init methods before calling this method.");
		}
	}

	/*
	 * Next, determine storage format and allocate space for solution set.
	 * For a more complete description of the arguments and what they are,
	 * how big they need to be, etc. please see the docs on DGBSV in LAPACK.
	 */
	__CLPK_integer		n = [self getRowCount] * [self getColCount];
	__CLPK_integer		kl = MIN([self getRowCount], [self getColCount]);
	__CLPK_integer		ku = kl;
	__CLPK_integer		klpku = kl + ku;
	BOOL				rowMajor = (kl == [self getColCount]);
	__CLPK_integer		nrhs = 1;
	__CLPK_integer		ldab = 2*kl + ku + 1;
	__CLPK_doublereal	*ab = NULL;
	if (!error) {
		ab = (__CLPK_doublereal	*) malloc( ldab*n*sizeof(__CLPK_doublereal) );
		if (ab == NULL) {
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - while trying to allocate the banded A matrix storage (%dx%d) for the solution, we ran into an allocation problem and couldn't get it. Please check into this as soon as possilbe.", ldab, n);
		}
	}
	__CLPK_integer		ldb = n;
	__CLPK_doublereal	*b = NULL;
	if (!error) {
		b = (__CLPK_doublereal *) malloc( ldb*nrhs*sizeof(__CLPK_doublereal) );
		if (b == NULL) {
			// free up what we've already obtained
			free(ab);
			// ...and log the error
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - while trying to allocate the RHS b matrix storage (%dx%d) for the solution, we ran into an allocation problem and couldn't get it. Please check into this as soon as possilbe.", ldb, nrhs);
		}
	}
	__CLPK_integer		*ipiv = NULL;
	if (!error) {
		ipiv = (__CLPK_integer *) malloc( n*sizeof(__CLPK_integer) );
		if (ipiv == NULL) {
			// free up what we've already obtained
			free(b);
			free(ab);
			// ...and log the error
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - while trying to allocate the ipivot storage (%dx1) for the solution, we ran into an allocation problem and couldn't get it. Please check into this as soon as possilbe.", n);
		}
	}
	
	/*
	 * Now we can populate the matricies with the equations to solve. We
	 * do this by going through all the nodes and placing the equation for
	 * that node into the cLAPACK matricies ab[] and b[] based on the
	 * location of each node, and the physical parameters for that node.
	 *
	 * Because cLAPACK uses banded storage and a vector as opposed to a
	 * two-dimensional array of numbers, we need to be able to convert
	 * the matrix notation into the actual values placed in the banded
	 * storage. Here's how it goes:
	 *
	 * The banded storage format is given as:
	 *     ab(kl+ku+1+i-j, j) = a(i, j)
	 * where i, j are constrained to be in the banded storage "coverage"
	 * and are both in the range (1,n), i.e. FORTRAN-style. To convert
	 * this to C-style arrays:
	 *     ab[kl+ku+i-j][j] = a[i][j]
	 * where now i and j run from 0..(n-1). The important difference is
	 * in the loss of the "+1" for the limiting case of kl=ku=0 needs to
	 * map to row=0 not row=1. The conversion to the row-major storage is,
	 * in general, accomplished with:
	 *     ab[col*ldab + row] = ab[row][col]
	 * for a given row and column, where again, row and col are in the range
	 * (0,n-1). Therefore, to map a general FORTRAN-style i,j from
	 * A into the banded storage is:
	 *     ab[j*ldab + kl+ku+i-j] = a(i, j)
	 * with the addition of a simple variable for speed:
	 *     ab[j*ldab + klpku + i-j] = a(i, j)
	 * and that's what we're implementing. B is implemented in a similar
	 * manner where:
	 *     b[rhs*ldb + i] = b(i, rhs)
	 * or:
	 *     b[i] = b(i)
	 * assuming that nrhs = 1, which it does for all our work.
	 */
	if (!error) {
		int		row = 0;
		int		col = 0;
		int		rows = [self getRowCount];
		int		cols = [self getColCount];
		// these will be the node numbers in the simulation
		int		ijn = 0;
		int		ln = 0;
		int		rn = 0;
		int		tn = 0;
		int		bn = 0;
		// these are the components of the coefficients
		__CLPK_doublereal	invHx2 = 1.0/([self getDeltaX] * [self getDeltaX]);
		__CLPK_doublereal	invHy2 = 1.0/([self getDeltaY] * [self getDeltaY]);
		// these are the values of rho and er at the node in the simulation
		__CLPK_doublereal	rho = 0;
		__CLPK_doublereal	er = 0;
		// now loop through all the node in the workspace
		for (row = 0; row < rows; row++) {
			for (col = 0; col < cols; col++) {
				/*
				 * We need to convert the row and col into a node number based
				 * on the most efficient banding possible.
				 */
				ijn = (rowMajor ? (row * cols + col) : (col * rows + row));
				
				// see if there's a fixed potential at this node
				if ([[self getVoltage] haveValueAtRow:row andCol:col]) {
					/*
					 * In this case the equation is simple: v(i,j) = v_fixed
					 */
					ab[ijn*ldab + klpku] = 1.0;
					b[ijn] = [self getVoltageAtNodeRow:row andCol:col];
				} else {
					/*
					 * In this case, we need to build up the entire Poission's Eq.
					 * for this node.
					 */
					// first, do the 'ij' node
					ab[ijn*ldab + klpku] = -2.0*(invHx2 + invHy2);
					// next, do the 'top' node
					if (row == 0) {
						// due to symmetry, add to the 'bottom' node
						bn = (rowMajor ? (ijn + cols) : (ijn + 1));
						ab[bn*ldab + klpku + ijn-bn] += invHy2;
					} else {
						tn = (rowMajor ? (ijn - cols) : (ijn - 1));
						ab[tn*ldab + klpku + ijn-tn] += invHy2;
					}
					// next, do the 'bottom' node
					if (row == (rows - 1)) {
						// due to symmetry, add to the 'top' node
						tn = (rowMajor ? (ijn - cols) : (ijn - 1));
						ab[tn*ldab + klpku + ijn-tn] += invHy2;
					} else {
						bn = (rowMajor ? (ijn + cols) : (ijn + 1));
						ab[bn*ldab + klpku + ijn-bn] += invHy2;
					}
					// next, do the 'left' node
					if (col == 0) {
						// due to symmetry, add to the 'right' node
						rn = (rowMajor ? (ijn + 1) : (ijn + rows));
						ab[rn*ldab + klpku + ijn-rn] += invHx2;
					} else {
						ln = (rowMajor ? (ijn - 1) : (ijn - rows));
						ab[ln*ldab + klpku + ijn-ln] += invHx2;
					}
					// finally, do the 'right' node
					if (col == (cols - 1)) {
						// due to symmetry, add to the 'left' node
						ln = (rowMajor ? (ijn - 1) : (ijn - rows));
						ab[ln*ldab + klpku + ijn-ln] += invHx2;
					} else {
						rn = (rowMajor ? (ijn + 1) : (ijn + rows));
						ab[rn*ldab + klpku + ijn-rn] += invHx2;
					}
					
					/*
					 * Now let's calculate the RHS of Ax=b...
					 */
					rho = [self getRhoAtNodeRow:row andCol:col];
					er = [self getEpsilonRAtNodeRow:row andCol:col];
					b[ijn] = -1.0 * rho/(er == 0 ? 1.0 : er);
				}
			}
		}
	}

	// finally, we can solve the system for the user using dgbsv_ in cLAPACK
	if (!error) {
		__CLPK_integer	info = 0;
		dgbsv_(&n, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &ldb, &info);
		if (info < 0) {
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - argument #%d had an illegal value to DGBSV in LAPACK. Please check into this.", -1*info);
		} else if (info > 0) {
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - diagonal #%d is zero indicating singularity which shouldn't happen.", info);
		}
	}

	// we need to create the resultant voltage matrix and populate it
	MaskedMatrix*		rv = nil;
	if (!error) {
		rv = [[[MaskedMatrix alloc] initWithRows:[self getRowCount] andCols:[self getColCount]] autorelease];
		if (rv == nil) {
			error = YES;
			NSLog(@"[SimWorkspace -simulateWorkspace] - the resultant voltage matrix for the simulation could not be created and this is a serious storage problem. The request was made for a %dx%d sized matrix, and that seems to be too much. Check into this.", [self getRowCount], [self getColCount]);
		} else {
			// now fill in all the values from the solution set
			int		row = 0;
			int		col = 0;
			int		ij = 0;
			for (row = 0; row < [self getRowCount]; row++) {
				for (col = 0; col < [self getColCount]; col++) {
					/*
					 * We need to convert the row and col into a node number based
					 * on the most efficient banding possible.
					 */
					ij = (rowMajor ? (row * [self getColCount] + col) : (col * [self getRowCount] + row));
					// now save it
					[rv setValue:b[ij] atRow:row andCol:col];
				}
			}
					
			// ...and don't forget to save it for the user
			[self _setResultantVoltage:rv];
		}
	}
	
	// in the end, we can release what it is that we don't need
	if (ipiv != NULL) {
		free(ipiv);
	}
	if (b != NULL) {
		free(b);
	}
	if (ab != NULL) {
		free(ab);
	}

	return !error;
}