This is one the function which computes the optimal score of a sequence alignment. The cost of the deletion is affine, it corresponds to a maximum of an entire row or column adding a constant amount per entry. In simple terms, an affine deletion penalty can be defined as having two costs: the cost for starting a deletion and an incremental cost for each additional, contiguous, deletion. Code which is constant and not relevant has been removed for clarity.

The following are the main reasons for automatic code generation in this example.

• Several different flavours of dynamic programming are required. Also on proteins and on DNA.
• Local/Global alignments.
• The mathematical operations performed in this function are quite simple, and hence very simple to optimize. The work with pointers to the arrays, instead of indices substantially increases the speed in most computers. This is generated automatically by the Maple generator.
• The program also converts the scoring matrix into integers and does all the computation with the largest possible integers.
• Finding the alignment will use two other similar versions of the code.

FUN(ALGEB ) LocalAlign ( Match *m, DayMatrix *dm )

{ iDayMatrix *idm;

  . . . . . .

  Mi1 = t1;
  Mi = t2;
  for( i=0; i <= lA; i++ ) {
	Mi1[i] = Mi[i] = 0;
	delcol[i] = -999999999;
	}
  imax = jmax = 1;  Mimaxjmax = 0;

  for( i=1; i <= lB; i++ ) {
    delrow = -999999999;
    Mij = 0;
    chB = idm->Sim[AToInt(s2+i-1)];
    delcolj = delcol;
    s2s1i = s2-s1+i;

    for( j=1; j <= lA; j++ ) {

	delcolj++;
	if( *delcolj >= Mi1[j]+diffDel )
	     *delcolj += idmDelIncr;
	else *delcolj = Mi1[j]+idmDelFixed;

	if( delrow >= Mij+diffDel )
	     delrow += idmDelIncr;
	else delrow = Mij+idmDelFixed;

	Mij = j==s2s1i ? -999999999 : chB[APs1[j-1]] + Mi1[j-1];

	if( *delcolj > delrow )
	     { if( *delcolj > Mij ) Mij = *delcolj; }
	else if( delrow > Mij ) Mij = delrow;

	if( Mij < 0 ) Mij = 0;
	else if( Mij > Mimaxjmax )
	    { Mimaxjmax = Mij;  imax = i;  jmax = j; }
	Mi[j] = Mij;
	}
    t3 = Mi1;  Mi1 = Mi;  Mi = t3;  /* swap Mi and Mi1 */
    }

  if( Mimaxjmax==0 ) goto nullreturn;

  lA = jmax;  lB = imax;
  for( i=0; i <= lA; i++ ) {
	Mi1[i] = Mi[i] = 0;
	delcol[i] = -999999999;
	}

  for( i=1; i <= lB; i++ ) {
    delrow = -999999999;
    Mij = 0;
    chB = idm->Sim[AToInt(s2+lB-i)];
    delcolj = delcol;
    s2s1i = s1-s2+lA-lB+i;

    for( j=1; j <= lA; j++ ) {

	delcolj++;
        if( *delcolj >= Mi1[j]+diffDel )
             *delcolj += idmDelIncr;
        else *delcolj = Mi1[j]+idmDelFixed;

        if( delrow >= Mij+diffDel )
             delrow += idmDelIncr;
        else delrow = Mij+idmDelFixed;

        Mij = j==s2s1i ? -999999999 : chB[APs1[lA-j]] + Mi1[j-1];

        if( *delcolj > delrow )
             { if( *delcolj > Mij ) Mij = *delcolj; }
        else if( delrow > Mij ) Mij = delrow;

        if( Mij >= Mimaxjmax ) {
	    if( Mij > Mimaxjmax ) userror("LocalAlign: should not happen");
	    m = (Match*)copyalg(A(m));
	    m->Sim = (double)Mimaxjmax / idm->ifac;
	    m->Offset1 = s1 - DB->string + lA-j;
	    m->Offset2 = s2 - DB->string + lB-i;
	    m->Length1 = j;
	    m->Length2 = i;
	    m->PamNumber = dm->PamNumber;
	    return( A(m) );
	    }
	Mi[j] = Mij;
        }
    t3 = Mi1;  Mi1 = Mi;  Mi = t3;  /* swap Mi and Mi1 */
    }
  userror("LocalAlign: Assertion failed");
  return(0);

}