diff options
Diffstat (limited to 'src/matrixOperations/division/sdiva.c')
| -rw-r--r-- | src/matrixOperations/division/sdiva.c | 167 |
1 files changed, 167 insertions, 0 deletions
diff --git a/src/matrixOperations/division/sdiva.c b/src/matrixOperations/division/sdiva.c new file mode 100644 index 00000000..a20d14df --- /dev/null +++ b/src/matrixOperations/division/sdiva.c @@ -0,0 +1,167 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-2008 - INRIA - Allan SIMON + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + + +/**** WARNING NOT WORK AT ALL FOR THE MOMENT ***/ +/**** Because of problem of conversion float-> double ****/ +#include "matrixDivision.h" +#include "lapack.h" +#include "stdio.h" + +void srdiva ( float* in1, int lines1, int columns1 , + float* in2, int lines2, int columns2 , + float* out ){ + +/* creation of all temp variable , maybe some of them are not accurate and could be deleted +**/ + int iexit = 0 ; + int iwork = 0 ; + int info = 0 ; + int imax = 0 ; + /*int rank = 0 ;*/ + + int tMinLinCol1 = 0 ; + int temp = 0; + + int i = 0 , j = 0 ,ij = 0 ,ji = 0; + + + int* pIpiv = (int*) malloc( sizeof(int) * (unsigned int) columns1); + int* pIwork = (int*) malloc( sizeof(int) * (unsigned int) columns1); + /*int* pJpvt = (int*) malloc( sizeof(int) * (unsigned int) lines1);*/ + + double* work = NULL; + double anorm = 0; + double rcond = 0 ; + double epsilon = getRelativeMachinePrecision() ; + + + + + float* transpOfIn1 = (float*) malloc(sizeof(float) * (unsigned int) lines1 * (unsigned int) columns1); + float* transpOfIn2 = (float*) malloc(sizeof(float) * (unsigned int) lines2 * (unsigned int) columns2); + + float* copyOfTransIn1 = (float*) malloc( sizeof(float) * (unsigned int) lines1 * (unsigned int) columns1); + + + +printf("\n>fin alloc \n"); + if (lines1 > columns1) + tMinLinCol1 = columns1 ; + else + tMinLinCol1 = lines1; + + if ( (tMinLinCol1 + 3 * lines1 + 1) > (2*tMinLinCol1 + lines2)) + temp = tMinLinCol1 + 3 * lines1 + 1; + else + temp = 2*tMinLinCol1 + lines2; + + if ( 4* columns1 > temp ) + iwork = 4* columns1 ; + else + iwork = temp ; + + work = (double*) malloc (sizeof(double) *(unsigned int) iwork ); + printf("\n>debut partie rigolote \n"); + anorm = getOneNorm( &lines1, &columns1 , (double* ) in1 , work); + printf("\n>debut partie rigolote 2\n"); + +/* end of allocation area*/ + + + stransposea ( in1, lines1, columns1, transpOfIn1); + stransposea ( in2, lines2, columns2, transpOfIn2) ; + + if ( lines1 == columns1 ) + { + srowcata ( in1, lines1, columns1, NULL, 0 , 0 , copyOfTransIn1 ) ; + + + /*/ put here algo of LU fact of in1 dgetrf ( &columns1 , &columns1 , in1 , &columns1 , pIpiv , &info ) + //return value in pIpiv*/ + dgetrf_ ( &columns1 , &columns1 , (double*) in1 , &columns1 , pIpiv , &info ); + + if ( info == 0 ) + { + /*/to get 1-norm of in1 put here algo dgecon which return rcond*/ + dgecon_("1" , &columns1, (double*) copyOfTransIn1 , &columns1, &anorm, + &rcond, work, pIwork, &info); + + if ( rcond > sqrt(epsilon )) + {/* put here algo to resolv linear equation in1 * X = in2 , the return value go in in2 + // put here algo to copy in2 in out */ + resolveSystemLinear (&columns1, &lines2,(double*) copyOfTransIn1, pIpiv, + (double* )transpOfIn2, &info) ; + stransposea ( transpOfIn2 , columns2 , lines2 , out) ; + iexit = 1; + } + + } + + + } + + if ( iexit== 0) + { + rcond = sqrt(epsilon ); + if (lines1 < columns1) + imax = columns1 ; + else + imax = lines1; + + +/* cette fonction genere des nombres aleatoires , youpi ! */ + /* dgelsy_ (&columns1, &lines1, &lines2, (double*) transpOfIn1 , &columns1, + (double*) transpOfIn2, &imax, pJpvt, &rcond, &rank, work, + &iwork, &info);*/ + + + + for(j = 0 ; j < lines1 ; j++) + { + for(i = 0 ; i < lines2 ; i++) + { + ij = i + j * lines2; + if ( lines1 > columns1 ) + ji = j + i *lines1 ; + else + ji = j + i *columns1; + out[ij] = transpOfIn2[ji]; + } + } + + + /* here put algo dgelsy1 which do maybe interesting things but i dunno what + + // if info return by dgelsy1 != 0 + // return ;*/ + + + + } + + printf("\n>debut partie 2\n"); + + + /* then we free all the allocated memory */ + + free ( pIpiv ) ; + free ( pIwork) ; + free ( work ) ; + + /* le probleme vient de ces trois free */ + free (transpOfIn1) ; + free (transpOfIn2) ; + free (copyOfTransIn1); + +} + |