; Performs the PCA decomposition of a dataset, returning the eigenvalues and the
; eigenvectors
; INPUT
;  A: NxM matrix that contains, as rows, the N M-dimensional observations
;    If the dataset are images of dimensions PxQ, we put each row as a M-dimensional vector
;    of size P*Q
; OUTPUT
;  eigenvalues: the set of M eigenvalues
;  eigenvectors: a matrix of size MxM containing the eigenvectors
;        eigenvectors[0,*] is the first eigenvector, associated with eigenvalues[0]
;        eigenvectors[i,*] is the i-th eigenvector, associated with eigenvalues[i]
;  
pro mypca, A_input, eigenvalues, eigenvectors, norm=norm, mincomp99=mincomp99, limite=limite, outmincomp=outmincomp
        A=A_input
        if keyword_set(norm) then begin
           pmess,'Using feature scaling: (x-mean)/stddev(x)'
           dim=size(A)
           s=fltarr(dim[1])
           for j=0, dim[1]-1 do begin
               aux=moment(A[j,*], mean=colmean, sdev=colsdev)
               A[j,*]=(A[j,*]-colmean)/colsdev
               ;colmean=total(A[*,j])/float(dim[1])
               ;a[*,j]=(a[*,j]-colmean)
               ;s[j]=sqrt(total(a[*,j]*a[*,j]))/float(dim[1])
               ;a[*,j]=a[*,j]/s[j]
           endfor
        endif

        ;if keyword_set(norm) then begin
        ;   pmess,'Using feature scaling: (x-mean)/stddev(x)'
        ;   dim=size(A)
        ;   s=fltarr(dim[2])
        ;   for j=0, dim[2]-1 do begin
        ;       aux=moment(A[*,j], mean=colmean, sdev=colsdev)
        ;       A[*,j]=(A[*,j]-colmean)/colsdev
        ;       ;colmean=total(A[*,j])/float(dim[1])
        ;       ;a[*,j]=(a[*,j]-colmean)
        ;       ;s[j]=sqrt(total(a[*,j]*a[*,j]))/float(dim[1])
        ;       ;a[*,j]=a[*,j]/s[j]
        ;   endfor
        ;endif
        ;;stop
	corr = matrix_multiply(A,A,/ATranspose)
	;corr = transpose(A)##A
	;svdc, corr, w, u, v
	la_svd, corr, w, u, v, status=converge
        if converge ne 0 then print,'PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS PROBLEMAS'  
	eigenvectors = v
	eigenvalues = w
        mincomp=fltarr((size(w))[1])
        for j=0,(size(w))[1]-1 do mincomp[j]=total(w[0:j])/total(w)
        outmincomp=mincomp
        lim0=0.99
        if keyword_set(limite) then lim0=limite
        mincomp99=min(where(mincomp ge lim0))
        mincomp99=[mincomp99, mincomp[mincomp99]]
end