Max 5 API Reference

jit.la.inverse.c

/* 
    Copyright 2001-2005 - Cycling '74
    Joshua Kit Clayton jkc@cycling74.com    
*/

#include "jit.common.h"

#define JIT_LA_INVERSE_DEF_THRESH 0.000000001

typedef struct _jit_la_inverse 
{
    t_object    ob;
    double      thresh;
} t_jit_la_inverse;

void *_jit_la_inverse_class;

t_jit_la_inverse *jit_la_inverse_new(void);
void jit_la_inverse_free(t_jit_la_inverse *x);
t_jit_err jit_la_inverse_matrix_calc(t_jit_la_inverse *x, void *inputs, void *outputs);
void jit_la_inverse_ident(t_jit_matrix_info *in_minfo, char *bip);

long jit_la_inverse_calculate_normupper_2d_real(t_jit_la_inverse *x, long *swapcount, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp);
long jit_la_inverse_calculate_normupper_2d_complex(t_jit_la_inverse *x, long *swapcount, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp);
void jit_la_inverse_calculate_backsub_2d_real(t_jit_la_inverse *x, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp);
void jit_la_inverse_calculate_backsub_2d_complex(t_jit_la_inverse *x, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp);

t_jit_err jit_la_inverse_init(void) 
{
    long attrflags=0;
    t_jit_object *attr,*mop,*o;
    t_atom a[2];
        
    _jit_la_inverse_class = jit_class_new("jit_la_inverse",(method)jit_la_inverse_new,(method)jit_la_inverse_free,
        sizeof(t_jit_la_inverse),0L);

    //add mop
    mop = jit_object_new(_jit_sym_jit_mop,1,1);
    o = jit_object_method(mop,_jit_sym_getoutput,1);
    jit_attr_setlong(o,_jit_sym_maxplanecount,2);
    jit_atom_setsym(a,_jit_sym_float32); //default
    jit_atom_setsym(a+1,_jit_sym_float64); 
    jit_object_method(o,_jit_sym_types,2,a);
    jit_class_addadornment(_jit_la_inverse_class,mop);
    //add methods
    jit_class_addmethod(_jit_la_inverse_class, (method)jit_la_inverse_matrix_calc,      "matrix_calc",      A_CANT, 0L);
    //add attributes    
    attrflags = JIT_ATTR_GET_DEFER_LOW | JIT_ATTR_SET_USURP_LOW;
    attr = jit_object_new(_jit_sym_jit_attr_offset,"thresh",_jit_sym_float64,attrflags,
        (method)0L,(method)0L,calcoffset(t_jit_la_inverse,thresh));
    jit_class_addattr(_jit_la_inverse_class,attr);
        
    jit_class_register(_jit_la_inverse_class);

    return JIT_ERR_NONE;
}

t_jit_err jit_la_inverse_matrix_calc(t_jit_la_inverse *x, void *inputs, void *outputs)
{
    t_jit_err err=JIT_ERR_NONE;
    long in_savelock,out_savelock;
    t_jit_matrix_info in_minfo,out_minfo,tmp_minfo;
    char *in_bp,*out_bp,*tmp_bp;
    long i,dimcount,planecount,dim[JIT_MATRIX_MAX_DIMCOUNT];
    long tmpsize,swapcount;
    void *tmp_matrix=NULL;
    void *in_matrix,*out_matrix;
    
    in_matrix   = jit_object_method(inputs,_jit_sym_getindex,0);
    out_matrix  = jit_object_method(outputs,_jit_sym_getindex,0);
    
    if (x&&in_matrix&&out_matrix) {
        swapcount=0;
        in_savelock = (long) jit_object_method(in_matrix,_jit_sym_lock,1);
        out_savelock = (long) jit_object_method(out_matrix,_jit_sym_lock,1);
        jit_object_method(in_matrix,_jit_sym_getinfo,&in_minfo);
        jit_object_method(out_matrix,_jit_sym_getinfo,&out_minfo);
        jit_object_method(in_matrix,_jit_sym_getdata,&in_bp);
        jit_object_method(out_matrix,_jit_sym_getdata,&out_bp);
        if (!in_bp) { err=JIT_ERR_INVALID_INPUT; goto out;}
        if (!out_bp) { err=JIT_ERR_INVALID_OUTPUT; goto out;}
        //compatible types?
        if (((in_minfo.type!=_jit_sym_float32)&&(in_minfo.type!=_jit_sym_float64))||
            (in_minfo.type!=out_minfo.type)) { err=JIT_ERR_MISMATCH_TYPE; goto out;}        
        //compatible planecounts?
        //planecount 1=real. planecount 2=complex
        if (((in_minfo.planecount!=1)&&(in_minfo.planecount!=2))||
            (in_minfo.planecount!=out_minfo.planecount)) { err=JIT_ERR_MISMATCH_PLANE; goto out;}       
        //check dimensions (square?)
        if ((in_minfo.dim[0]!=out_minfo.dim[1])||
            (in_minfo.dim[0]!=out_minfo.dim[0])||
            (in_minfo.dim[1]!=out_minfo.dim[1])) { err=JIT_ERR_MISMATCH_DIM; goto out;}
        //get dimensions/planecount
        dimcount   = out_minfo.dimcount;
        planecount = out_minfo.planecount;          
        
        //set output matrix to be identity
        jit_object_method(out_matrix,_jit_sym_clear);
        jit_la_inverse_ident(&out_minfo, out_bp);
        
        //first copy matrix, and then calculate normalized upper trianglular matrix inplace 
        if (!(tmp_matrix=jit_object_new(_jit_sym_jit_matrix,&in_minfo)))
            { err=JIT_ERR_OUT_OF_MEM; goto out;}
        jit_object_method(tmp_matrix,_jit_sym_lock,1);
        jit_object_method(tmp_matrix,_jit_sym_frommatrix,in_matrix,NULL);
        jit_object_method(tmp_matrix,_jit_sym_getinfo,&tmp_minfo);
        jit_object_method(tmp_matrix,_jit_sym_getdata,&tmp_bp);
        if (!tmp_bp) goto cleanup;
        
        if (planecount==1) {
            if (jit_la_inverse_calculate_normupper_2d_real(x, &swapcount, &tmp_minfo, tmp_bp, &out_minfo, out_bp))
                { jit_object_method(out_matrix,_jit_sym_clear); err = JIT_ERR_GENERIC; goto cleanup; }
            jit_la_inverse_calculate_backsub_2d_real(x, &tmp_minfo, tmp_bp, &out_minfo, out_bp);
        } else {
            if (jit_la_inverse_calculate_normupper_2d_complex(x, &swapcount, &tmp_minfo, tmp_bp, &out_minfo, out_bp))
                { jit_object_method(out_matrix,_jit_sym_clear); err = JIT_ERR_GENERIC; goto cleanup; }  
            jit_la_inverse_calculate_backsub_2d_complex(x, &tmp_minfo, tmp_bp, &out_minfo, out_bp);
        }
        
cleanup:
        if (tmp_matrix) {
            jit_object_method(tmp_matrix,_jit_sym_lock,0);
            jit_object_free(tmp_matrix);
        }
    } else {
        return JIT_ERR_INVALID_PTR;
    }
    
out:
    jit_object_method(out_matrix,_jit_sym_lock,out_savelock);
    jit_object_method(in_matrix,_jit_sym_lock,in_savelock);
    return err;
}

//assumes input is all zeros
void jit_la_inverse_ident(t_jit_matrix_info *in_minfo, char *bip)
{
    long i,rows;
    char *p;

    rows = in_minfo->dim[1];    
    if (in_minfo->type==_jit_sym_float32) {
        for (i=0;i<rows;i++) {
            p = bip  + i*in_minfo->dimstride[1] + i*in_minfo->dimstride[0]; 
            *((float *)p) = 1.;
        }
    } else if (in_minfo->type==_jit_sym_float64) {
        for (i=0;i<rows;i++) {
            p = bip  + i*in_minfo->dimstride[1] + i*in_minfo->dimstride[0]; 
            *((double *)p) = 1.;
        }
    } 
}


long jit_la_inverse_calculate_normupper_2d_real(t_jit_la_inverse *x, long *swapcount, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp)
{
    long i,j,n,rows;
    t_jit_op_info in_opinfo,tmp_opinfo,p_opinfo,p_tmp_opinfo;
    double val,val2,val_inv;
    char *op;

    n = in_minfo->dim[0];   
    if (in_minfo->type==_jit_sym_float32) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/4; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/4;
        rows = in_minfo->dim[1];
        for (i=0;i<rows;i++) {
            in_opinfo.p = bip  + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            val = *((float *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]));
            //perform row swapping if necessary             
            if (ABS(val)<x->thresh) { 
                j=i+1;
                while ((j<rows)&&(ABS(val)<x->thresh)) { 
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1];     
                    val = *((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                    j++;
                }
                if (!(ABS(val)<x->thresh)) {
                    jit_op_vector_swap_float32(n,NULL,&in_opinfo,&tmp_opinfo,NULL);
                    jit_op_vector_swap_float32(n,NULL,&p_opinfo,&p_tmp_opinfo,NULL);
                    *swapcount += 1;
                } else {
                    return -1;
                }
            }
            
            if (!(ABS(val)<x->thresh)) {  // thresh?
                //normalize
                val2 = 1./val;
                jit_op_vector_ax_float32(n,&val2,&in_opinfo,NULL,&in_opinfo); // note in place
                jit_op_vector_ax_float32(n,&val2,&p_opinfo,NULL,&p_opinfo); // note in place                    
                //eliminate all non zero values in column i
                j=i+1;
                val_inv=-1.;
                while (j<rows) {
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                    val2 = *((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                    if (!(ABS(val2)<x->thresh)) { // thresh?
                        val2 = val2*val_inv;
                        jit_op_vector_axpy_float32(n,&val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                        jit_op_vector_axpy_float32(n,&val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                    }
                    j++;
                }
            } 
        }
    } else if (in_minfo->type==_jit_sym_float64) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/8; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/8;
        rows = in_minfo->dim[1];
        for (i=0;i<rows;i++) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            val = *((double *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]));
            //perform row swapping if necessary             
            if ((ABS(val)<x->thresh)) { 
                j=i+1;
                while ((j<rows)&&(ABS(val)<x->thresh)) { 
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1];     
                    val = *((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                    j++;
                }
                if (!(ABS(val)<x->thresh)) {
                    jit_op_vector_swap_float64(n,NULL,&in_opinfo,&tmp_opinfo,NULL);
                    jit_op_vector_swap_float64(n,NULL,&p_opinfo,&p_tmp_opinfo,NULL);
                    *swapcount += 1;
                } else {
                    return -1;
                }
            }
            if (!(ABS(val)<x->thresh)) {  // thresh?
                //normalize
                val2 = 1./val;
                jit_op_vector_ax_float64(n,&val2,&in_opinfo,NULL,&in_opinfo); // note in place
                jit_op_vector_ax_float64(n,&val2,&p_opinfo,NULL,&p_opinfo); // note in place                    
                //eliminate all non zero values in column i
                j=i+1;
                val_inv=-1.;
                while (j<rows) {
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                    val2 = *((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                    if (!(ABS(val2)<x->thresh)) { // thresh?
                        val2 = val2*val_inv;
                        jit_op_vector_axpy_float64(n,&val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                        jit_op_vector_axpy_float64(n,&val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                    }
                    j++;
                }
            }
        }
    }
    return 0; 
}

void jit_la_inverse_calculate_backsub_2d_real(t_jit_la_inverse *x, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp)
{
    long i,j,n,rows;
    t_jit_op_info in_opinfo,tmp_opinfo,p_opinfo,p_tmp_opinfo;
    double val2;
    char *op;

    n = in_minfo->dim[0];   
    if (in_minfo->type==_jit_sym_float32) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/4; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/4;
        rows = in_minfo->dim[1];
        for (i=rows-1;i>=0;i--) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            
            //eliminate all non zero values in column i
            j=i-1;
            while (j>=0) {
                tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                val2 = *((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                if (!(ABS(val2)<x->thresh)) { // thresh?
                    val2 = -val2;
                    jit_op_vector_axpy_float32(n,&val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                    jit_op_vector_axpy_float32(n,&val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                }
                j--;
            }
        }
    } else if (in_minfo->type==_jit_sym_float64) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/8; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/8;
        rows = in_minfo->dim[1];
        for (i=rows-1;i>=0;i--) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            
            //eliminate all non zero values in column i
            j=i-1;
            while (j>=0) {
                tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                val2 = *((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]));
                if (!(ABS(val2)<x->thresh)) { // thresh?
                    val2 = -val2;
                    jit_op_vector_axpy_float64(n,&val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                    jit_op_vector_axpy_float64(n,&val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                }
                j--;
            }
        }
    } 
}

#define COMPLEX_IS_ZERO(a,b) (((a)==0)&&((b)==0))
#define COMPLEX_IS_ZERO_THRESH(a,b,thresh) ((ABS(a)<(thresh))&&(ABS(b)<(thresh)))
#define COMPLEX_RECIP_REAL(a,b) ((a)/(((a)*(a))+((b)*(b))))
#define COMPLEX_RECIP_IMAG(a,b) (-(b)/(((a)*(a))+((b)*(b))))
#define COMPLEX_MULT_REAL(a,b,c,d) (((a)*(c))-((b)*(d)))
#define COMPLEX_MULT_IMAG(a,b,c,d) (((a)*(d))+((b)*(c)))

long jit_la_inverse_calculate_normupper_2d_complex(t_jit_la_inverse *x, long *swapcount, t_jit_matrix_info *in_minfo, char *bip,t_jit_matrix_info *p_minfo, char *bpp )
{
    long i,j,n,rows;
    t_jit_op_info in_opinfo,tmp_opinfo,p_opinfo,p_tmp_opinfo;
    double val[2],val2[2],val3[2],val_inv[2];
    char *op;

    n = in_minfo->dim[0];   
    if (in_minfo->type==_jit_sym_float32) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/4; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/4;
        rows = in_minfo->dim[1];
        for (i=0;i<rows;i++) {
            in_opinfo.p = bip  + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            val[0] = ((float *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]))[0];
            val[1] = ((float *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]))[1];
            
            //perform row swapping if necessary             
            if (COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) { 
                j=i+1;
                while ((j<rows)&&COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) { 
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1];     
                    val[0] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                    val[1] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                    j++;
                }
                if (!COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) {
                    jit_op_vector_swap_float32_complex(n,NULL,&in_opinfo,&tmp_opinfo,NULL);
                    jit_op_vector_swap_float32_complex(n,NULL,&p_opinfo,&p_tmp_opinfo,NULL);
                    *swapcount += 1;
                } else {
                    return -1;
                }
            }
            if (!COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) {  // thresh?
                //normalize
                val2[0] = COMPLEX_RECIP_REAL(val[0],val[1]); 
                val2[1] = COMPLEX_RECIP_IMAG(val[0],val[1]);
                jit_op_vector_ax_float32_complex(n,val2,&in_opinfo,NULL,&in_opinfo); // note in place
                jit_op_vector_ax_float32_complex(n,val2,&p_opinfo,NULL,&p_opinfo); // note in place                 
                //eliminate all non zero values in column i
                j=i+1;
                val_inv[0] = -1; 
                val_inv[1] = -0;
                while (j<rows) {
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                    val2[0] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                    val2[1] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                    if (!COMPLEX_IS_ZERO_THRESH(val2[0],val2[1],x->thresh)) { // thresh?
                        val3[0] = COMPLEX_MULT_REAL(val2[0],val2[1],val_inv[0],val_inv[1]);
                        val3[1] = COMPLEX_MULT_IMAG(val2[0],val2[1],val_inv[0],val_inv[1]);
                        jit_op_vector_axpy_float32_complex(n,val3,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                        jit_op_vector_axpy_float32_complex(n,val3,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                    }
                    j++;
                }
            }
        }
    } else if (in_minfo->type==_jit_sym_float64) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/8; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/8;
        rows = in_minfo->dim[1];
        for (i=0;i<rows;i++) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            val[0] = ((double *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]))[0];
            val[1] = ((double *)(((char *)in_opinfo.p) + i*in_minfo->dimstride[0]))[1];
            
            //perform row swapping if necessary             
            if (COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) { 
                j=i+1;
                while ((j<rows)&&COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) { 
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1];     
                    val[0] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                    val[1] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                    j++;
                }
                if (!COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) {
                    jit_op_vector_swap_float64_complex(n,NULL,&in_opinfo,&tmp_opinfo,NULL);
                    jit_op_vector_swap_float64_complex(n,NULL,&p_opinfo,&p_tmp_opinfo,NULL);
                    *swapcount += 1;
                } else {
                    return -1;
                }
            }
            if (!COMPLEX_IS_ZERO_THRESH(val[0],val[1],x->thresh)) {  // thresh?
                //normalize
                val2[0] = COMPLEX_RECIP_REAL(val[0],val[1]); 
                val2[1] = COMPLEX_RECIP_IMAG(val[0],val[1]);
                jit_op_vector_ax_float64_complex(n,val2,&in_opinfo,NULL,&in_opinfo); // note in place
                jit_op_vector_ax_float64_complex(n,val2,&p_opinfo,NULL,&p_opinfo); // note in place                 
                //eliminate all non zero values in column i
                j=i+1;
                val_inv[0] = -1; 
                val_inv[1] = -0;
                while (j<rows) {
                    tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                    p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                    val2[0] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                    val2[1] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                    if (!COMPLEX_IS_ZERO_THRESH(val2[0],val2[1],x->thresh)) { // thresh?
                        val3[0] = COMPLEX_MULT_REAL(val2[0],val2[1],val_inv[0],val_inv[1]);
                        val3[1] = COMPLEX_MULT_IMAG(val2[0],val2[1],val_inv[0],val_inv[1]);
                        jit_op_vector_axpy_float64_complex(n,val3,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                        jit_op_vector_axpy_float64_complex(n,val3,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                    }
                    j++;
                }
            }
        }
    } 
    return 0;
}

void jit_la_inverse_calculate_backsub_2d_complex(t_jit_la_inverse *x, t_jit_matrix_info *in_minfo, char *bip, t_jit_matrix_info *p_minfo, char *bpp)
{
    long i,j,n,rows;
    t_jit_op_info in_opinfo,tmp_opinfo,p_opinfo,p_tmp_opinfo;
    double val2[2];
    char *op;

    n = in_minfo->dim[0];   
    if (in_minfo->type==_jit_sym_float32) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/4; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/4;
        rows = in_minfo->dim[1];
        for (i=rows-1;i>=0;i--) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            
            //eliminate all non zero values in column i
            j=i-1;
            while (j>=0) {
                tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                val2[0] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                val2[1] = ((float *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                if (!COMPLEX_IS_ZERO_THRESH(val2[0],val2[1],x->thresh)) { // thresh?
                    val2[0] = -val2[0];
                    val2[1] = -val2[1];
                    jit_op_vector_axpy_float32_complex(n,val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                    jit_op_vector_axpy_float32_complex(n,val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                }
                j--;
            }
        }
    } else if (in_minfo->type==_jit_sym_float64) {
        tmp_opinfo.stride = in_opinfo.stride = in_minfo->dimstride[0]/8; //row vector
        p_opinfo.stride = p_tmp_opinfo.stride = p_minfo->dimstride[0]/8;
        rows = in_minfo->dim[1];
        for (i=rows-1;i>=0;i--) {
            in_opinfo.p = bip + i*in_minfo->dimstride[1];
            p_opinfo.p = bpp + i*p_minfo->dimstride[1];     
            
            //eliminate all non zero values in column i
            j=i-1;
            while (j>=0) {
                tmp_opinfo.p = bip + j*in_minfo->dimstride[1]; //row j
                p_tmp_opinfo.p = bpp + j*p_minfo->dimstride[1]; 
                val2[0] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[0];
                val2[1] = ((double *)(((char *)tmp_opinfo.p) + i*in_minfo->dimstride[0]))[1];
                if (!COMPLEX_IS_ZERO_THRESH(val2[0],val2[1],x->thresh)) { // thresh?
                    val2[0] = -val2[0];
                    val2[1] = -val2[1];
                    jit_op_vector_axpy_float64_complex(n,val2,&in_opinfo,&tmp_opinfo,&tmp_opinfo); // note in place
                    jit_op_vector_axpy_float64_complex(n,val2,&p_opinfo,&p_tmp_opinfo,&p_tmp_opinfo); // note in place
                }
                j--;
            }
        }
    } 
}

t_jit_la_inverse *jit_la_inverse_new(void)
{
    t_jit_la_inverse *x;
        
    if (x=(t_jit_la_inverse *)jit_object_alloc(_jit_la_inverse_class)) {
        x->thresh = JIT_LA_INVERSE_DEF_THRESH; 
    } else {
        x = NULL;
    }   
    return x;
}

void jit_la_inverse_free(t_jit_la_inverse *x)
{
}

Copyright © 2008, Cycling '74