C = &alpha;A + &beta;u*v<sup>T</sup>

@param alpha scale factor applied to A
@param A matrix
@param beta scale factor applies to outer product
@param u vector
@param v vector
@param C Storage for solution. Can be same instance as A.
504830
public static void multAddOuter( double alpha , DMatrix5x5 A , double beta , DMatrix5 u , DMatrix5 v , DMatrix5x5 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a15 = alpha*A.a15 + beta*u.a1*v.a5;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a25 = alpha*A.a25 + beta*u.a2*v.a5;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a35 = alpha*A.a35 + beta*u.a3*v.a5;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
        C.a45 = alpha*A.a45 + beta*u.a4*v.a5;
        C.a51 = alpha*A.a51 + beta*u.a5*v.a1;
        C.a52 = alpha*A.a52 + beta*u.a5*v.a2;
        C.a53 = alpha*A.a53 + beta*u.a5*v.a3;
        C.a54 = alpha*A.a54 + beta*u.a5*v.a4;
        C.a55 = alpha*A.a55 + beta*u.a5*v.a5;
    }
1000736
public static void divide( DMatrix4 a , double alpha ) {
        a.a1 /= alpha;
        a.a2 /= alpha;
        a.a3 /= alpha;
        a.a4 /= alpha;
    }
1134771
def smeftpar(scale, C, basis):
    """"""
    # start with a zero dict and update it with the input values
    MW = p['m_W']
    # MZ = p['m_Z']
    GF = p['GF']
    Mh = p['m_h']
    vb = sqrt(1 / sqrt(2) / GF)
    v = vb  # TODO
    _d = vMh2_to_m2Lambda(v=v, Mh2=Mh**2, C=C)
    m2 = _d['m2'].real
    Lambda = _d['Lambda'].real
    gsbar = sqrt(4 * pi * p['alpha_s'])
    gs = (1 - C['phiG'] * (v**2)) * gsbar
    gbar = 2 * MW / v
    g = gbar * (1 - C['phiW'] * (v**2))
    ebar = sqrt(4 * pi * p['alpha_e'])
    gp = get_gpbar(ebar, gbar, v, C)
    c = {}
    c['m2'] = m2
    c['Lambda'] = Lambda
    c['g'] = g
    c['gp'] = gp
    c['gs'] = gs
    K = ckmutil.ckm.ckm_tree(p['Vus'], p['Vub'], p['Vcb'], p['delta'])
    if basis == 'Warsaw':
        Mu = K.conj().T @ np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = np.diag([p['m_d'], p['m_s'], p['m_b']])
    elif basis == 'Warsaw up':
        Mu = np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = K @ np.diag([p['m_d'], p['m_s'], p['m_b']])
    else:
        raise ValueError("Basis '{}' not supported".format(basis))
    Me = np.diag([p['m_e'], p['m_mu'], p['m_tau']])
    c['Gd'] = Md / (v / sqrt(2)) + C['dphi'] * (v**2) / 2
    c['Gu'] = Mu / (v / sqrt(2)) + C['uphi'] * (v**2) / 2
    c['Ge'] = Me / (v / sqrt(2)) + C['ephi'] * (v**2) / 2
    return c
757325
public static FloatMatrix gemv(float alpha, FloatMatrix a,
			FloatMatrix x, float beta, FloatMatrix y) {
		if (false) {
			NativeBlas.sgemv('N', a.rows, a.columns, alpha, a.data, 0, a.rows, x.data, 0,
					1, beta, y.data, 0, 1);
		} else {
			if (beta != 0.0f) {
				for (int i = 0; i < y.length; i++)
					y.data[i] = beta * y.data[i];
			} else {
				for (int i = 0; i < y.length; i++)
					y.data[i] = 0.0f;
			}

	
			for (int j = 0; j < a.columns; j++) {
				float xj = x.get(j);
				if (xj != 0.0f) {
					for (int i = 0; i < a.rows; i++)
							y.data[i] += alpha * a.get(i, j) * xj;
				}
			}
		}
		return y;
	}
879956
def m_c(mcmc, scale, f, alphasMZ=0.1185, loop=3):
    """
    if scale == mcmc:
        return mcmc  # nothing to do
    _sane(scale, f)
    crd = rundec.CRunDec()
    alphas_mc = alpha_s(mcmc, 4, alphasMZ=alphasMZ, loop=loop)
    if f == 4:
        alphas_scale = alpha_s(scale, f, alphasMZ=alphasMZ, loop=loop)
        return crd.mMS2mMS(mcmc, alphas_mc, alphas_scale, f, loop)
    elif f == 3:
        crd.nfMmu.Mth = 1.3
        crd.nfMmu.muth = 1.3
        crd.nfMmu.nf = 4
        return crd.mH2mL(mcmc, alphas_mc, mcmc, crd.nfMmu, scale, loop)
    elif f == 5:
        crd.nfMmu.Mth = 4.8
        crd.nfMmu.muth = 4.8
        crd.nfMmu.nf = 5
        return crd.mL2mH(mcmc, alphas_mc, mcmc, crd.nfMmu, scale, loop)
    else:
        raise ValueError("Invalid input: f={}, scale={}".format(f, scale))
721457
public void factorize(boolean checkSymmetry) throws Exception{
		if (checkSymmetry && !Property.TWELVE.isSymmetric(Q)) {
			throw new Exception("Matrix is not symmetric");
		}
		
		if(this.rescaler != null){
			double[] cn_00_original = null;
			double[] cn_2_original = null;
			double[] cn_00_scaled = null;
			double[] cn_2_scaled = null;
			if(log.isDebugEnabled()){
				cn_00_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_original Q before scaling: " + ArrayUtils.toString(cn_00_original));
				cn_2_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_original Q before scaling : " + ArrayUtils.toString(cn_2_original));
			}
			//scaling the Q matrix, we have:
			//Q1 = U.Q.U[T] = U.L.L[T].U[T] = (U.L).(U.L)[T] 
			//and because U is diagonal it preserves the triangular form of U.L, so
			//Q1 = U.Q.U[T] = L1.L1[T] is the new Cholesky decomposition  
			DoubleMatrix1D Uv = rescaler.getMatrixScalingFactorsSymm(Q);
			if(log.isDebugEnabled()){
				boolean checkOK = rescaler.checkScaling(ColtUtils.fillSubdiagonalSymmetricMatrix(Q), Uv, Uv);
				if(!checkOK){
					log.warn("Scaling failed (checkScaling = false)");
				}
			}
			this.U = Uv;
			this.Q = ColtUtils.diagonalMatrixMult(Uv, Q, Uv);
			if(log.isDebugEnabled()){
				cn_00_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_scaled Q after scaling : " + ArrayUtils.toString(cn_00_scaled));
				cn_2_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_scaled Q after scaling  : " + ArrayUtils.toString(cn_2_scaled));
				
				if(cn_00_original[0] < cn_00_scaled[0] || cn_2_original[0] < cn_2_scaled[0]){
					//log.info("Q: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()));
					log.warn("Problematic scaling");
					//throw new RuntimeException("Scaling failed");
				}
			}
		}
		
		//double threshold = Math.pow(Utils.getDoubleMachineEpsilon(), 2);
		double threshold = Utils.getDoubleMachineEpsilon();
		this.LData = new double[dim][];

		for (int i = 0; i < dim; i++) {
			LData[i] = new double[i+1];
			//DoubleMatrix1D QI = Q.viewRow(i);
			double[] LDataI = LData[i];
			//j < i
			for (int j = 0; j < i; j++) {
				double[] LDataJ = LData[j];
				double sum = 0.0;
				for (int k = 0; k < j; k++) {
					sum += LDataI[k] * LDataJ[k];
				}
				LDataI[j] = 1.0 / LDataJ[j] * (Q.getQuick(i, j) - sum);
			}
			//j==i
			double sum = 0.0;
			for (int k = 0; k < i; k++) {
				sum += Math.pow(LDataI[k], 2);
			}
			double d = Q.getQuick(i, i) - sum;
			if(!(d > threshold)){
				throw new Exception("not positive definite matrix");
			}
			LDataI[i] = Math.sqrt(d);
		}
	}
1175078
def rcauchy(alpha, beta, size=None):
    

    return alpha + beta * np.tan(pi * random_number(size) - pi / 2.0)
1260618
public SparseVector scale(double alpha) {
		SparseVector a = this;
		SparseVector c = new SparseVector(N);
		for (int i : a.symbolTable) c.put(i, alpha * a.get(i));
		return c;
	}
73307
public static void addCol(Matrix A, int j, double c)
    {
        addCol(A, j, 0, A.rows(), c);
    }
1179920
public static void multTransA_vecCol(final int blockLength ,
                                         DSubmatrixD1 A , DSubmatrixD1 B ,
                                         DSubmatrixD1 C )
    {
        int widthA = A.col1 - A.col0;
        if( widthA > blockLength )
            throw new IllegalArgumentException("A is expected to be at most one block wide.");

        for( int j = B.col0; j < B.col1; j += blockLength ) {
            int widthB = Math.min( blockLength , B.col1 - j );

            int indexC = C.row0*C.original.numCols + (j-B.col0+C.col0)*widthA;

            for( int k = A.row0; k < A.row1; k += blockLength ) {
                int heightA = Math.min( blockLength , A.row1 - k );

                int indexA = k*A.original.numCols + A.col0*heightA;
                int indexB = (k-A.row0+B.row0)*B.original.numCols + j*heightA;

                if( k == A.row0 )
                    multTransABlockSet_lowerTriag(A.original.data,B.original.data,C.original.data,
                            indexA,indexB,indexC,heightA,widthA,widthB);
                else
                    InnerMultiplication_DDRB.blockMultPlusTransA(A.original.data,B.original.data,C.original.data,
                            indexA,indexB,indexC,heightA,widthA,widthB);
            }
        }
    }
227785
public static boolean cholU( DMatrix5x5 A ) {

        A.a11 = Math.sqrt(A.a11);
        A.a21 = 0;
        A.a31 = 0;
        A.a41 = 0;
        A.a51 = 0;
        A.a12 = (A.a12)/A.a11;
        A.a22 = Math.sqrt(A.a22-A.a12*A.a12);
        A.a32 = 0;
        A.a42 = 0;
        A.a52 = 0;
        A.a13 = (A.a13)/A.a11;
        A.a23 = (A.a23-A.a12*A.a13)/A.a22;
        A.a33 = Math.sqrt(A.a33-A.a13*A.a13-A.a23*A.a23);
        A.a43 = 0;
        A.a53 = 0;
        A.a14 = (A.a14)/A.a11;
        A.a24 = (A.a24-A.a12*A.a14)/A.a22;
        A.a34 = (A.a34-A.a13*A.a14-A.a23*A.a24)/A.a33;
        A.a44 = Math.sqrt(A.a44-A.a14*A.a14-A.a24*A.a24-A.a34*A.a34);
        A.a54 = 0;
        A.a15 = (A.a15)/A.a11;
        A.a25 = (A.a25-A.a12*A.a15)/A.a22;
        A.a35 = (A.a35-A.a13*A.a15-A.a23*A.a25)/A.a33;
        A.a45 = (A.a45-A.a14*A.a15-A.a24*A.a25-A.a34*A.a35)/A.a44;
        A.a55 = Math.sqrt(A.a55-A.a15*A.a15-A.a25*A.a25-A.a35*A.a35-A.a45*A.a45);
        return !UtilEjml.isUncountable(A.a55);
    }
1262574
public void scaleYResolution(double factor) {
        this.y_resolution = Math.max(1.0, this.y_resolution * factor);
        updateYResolution();
        updateUpperYValue();
        updateCanvas(true);
    }
472375
public static void scale( double alpha , DMatrix6x6 a ) {
        a.a11 *= alpha; a.a12 *= alpha; a.a13 *= alpha; a.a14 *= alpha; a.a15 *= alpha; a.a16 *= alpha;
        a.a21 *= alpha; a.a22 *= alpha; a.a23 *= alpha; a.a24 *= alpha; a.a25 *= alpha; a.a26 *= alpha;
        a.a31 *= alpha; a.a32 *= alpha; a.a33 *= alpha; a.a34 *= alpha; a.a35 *= alpha; a.a36 *= alpha;
        a.a41 *= alpha; a.a42 *= alpha; a.a43 *= alpha; a.a44 *= alpha; a.a45 *= alpha; a.a46 *= alpha;
        a.a51 *= alpha; a.a52 *= alpha; a.a53 *= alpha; a.a54 *= alpha; a.a55 *= alpha; a.a56 *= alpha;
        a.a61 *= alpha; a.a62 *= alpha; a.a63 *= alpha; a.a64 *= alpha; a.a65 *= alpha; a.a66 *= alpha;
    }
523198
public static void jacobian_Control4( DMatrixRMaj L_full ,
										  double beta[] , DMatrixRMaj A )
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2]; double b3 = beta[3];

		final double ld[] = L_full.data;

		for( int i = 0; i < 6; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2 +   ld[li+3]*b3;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+4]*b1 +   ld[li+5]*b2 +   ld[li+6]*b3;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+5]*b1 + 2*ld[li+7]*b2 +   ld[li+8]*b3;
			A.data[indexA++] =   ld[li+3]*b0 +   ld[li+6]*b1 +   ld[li+8]*b2 + 2*ld[li+9]*b3;
		}
	}
897027
def syr2(X, u, v, alpha = 1.0, beta = 1.0, reordered=False):
    
    """
    assert X.is_factor is False, "cspmatrix factor object"
    symb = X.symb
    n = symb.n
    snptr = symb.snptr
    snode = symb.snode
    blkval = X.blkval
    blkptr = symb.blkptr
    relptr = symb.relptr
    snrowidx = symb.snrowidx
    sncolptr = symb.sncolptr

    if symb.p is not None and reordered is False:
        up = u[symb.p]
        vp = v[symb.p]
    else:
        up = u
        vp = v
            
    for k in range(symb.Nsn):
        nn = snptr[k+1]-snptr[k]     
        na = relptr[k+1]-relptr[k] 
        nj = na + nn

        for i in range(nn): blas.scal(beta, blkval, n = nj-i, offset = blkptr[k]+(nj+1)*i)

        uk = up[snrowidx[sncolptr[k]:sncolptr[k+1]]]
        vk = vp[snrowidx[sncolptr[k]:sncolptr[k+1]]]

        blas.syr2(uk, vk, blkval, n = nn, offsetA = blkptr[k], ldA = nj, alpha = alpha)
        blas.ger(uk, vk, blkval, m = na, n = nn, offsetx = nn, offsetA = blkptr[k]+nn, ldA = nj, alpha = alpha)
        blas.ger(vk, uk, blkval, m = na, n = nn, offsetx = nn, offsetA = blkptr[k]+nn, ldA = nj, alpha = alpha)
        
    return
481978
@Override
    public Matrix multAdd(double alpha, Matrix B, Matrix C) {
        checkMultAdd(B, C);
        if (alpha == 0)
            return C;
        for (int i = 0; i < numRows; i++) {
            Node row = links.startOfRow(i);
            if (row != null)
                for (int j = 0; j < B.numColumns(); j++) {
                    Node node = row;
                    double v = 0;
                    while (node != null && node.row == i) {
                        v += (B.get(node.col, j) * node.val);
                        node = node.rowTail;
                    }
                    if (v != 0)
                        C.add(i, j, alpha * v);
                }
        }
        return C;
    }
646640
public static void fill( DMatrix3x3 a , double v  ) {
        a.a11 = v; a.a12 = v; a.a13 = v;
        a.a21 = v; a.a22 = v; a.a23 = v;
        a.a31 = v; a.a32 = v; a.a33 = v;
    }